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FORMATION OF Z - 6 QUASARS FROM HIERARCHICAL GALAXY MERGERS 
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ABSTRACT 

The discovery of luminous quasars at redshift z ~ 6 indicates the presence of supermassive black holes 
(SMBHs) of mass ~ 10 9 M Q when the Universe was less than one billion years old. This finding presents 
several challenges for theoretical models, because whether such massive objects can form so early in the A- 
cold dark matter (ACDM) cosmology, the leading theory for cosmic structure formation, is an open question. 
Furthermore, whether the formation process requires exotic physics such as super-Eddington accretion remains 
undecided. Here, we present the first multi-scale simulations that, together with a self-regulated model for 
the SMBH growth, produce a luminous quasar at z ~ 6.5 in the ACDM paradigm. We follow the hierarchical 
assembly history of the most massive halo in a ~ 3Gpc 3 volume, and find that this halo of ~ 8 x 10 12 M Q 
forming at z ~ 6.5 after several major mergers is able to reproduce a n umber of observed properties of SDSS 
J 1 148+525 1 , the most distant quasar detected at z = 6.42 dFan et al.l2003l) . Moreover, the SMBHs grow through 
gas accretion below the Eddington limit in a self -regulated manner owing to feedback. We find that the pro- 
genitors experience vigorous star formation (up to 10 4 M© yr" 1 ) preceding the major quasar phase such that the 
stellar mass of the quasar host reaches 10 I2 M Q at z ~ 6.5, consistent with observations of significant metal en- 
richment in SDSS J 1 148+525 1 . The merger remnant thus obeys similar MBH - Mb u ig e scaling relation observed 
locally as a consequence of coeval growth and evolution of the SMBH and its host galaxy. Our results provide 
a viable formation mechanism for z ~ 6 quasars in the standard ACDM cosmology, and demonstrate a com- 
mon, merger-driven origin for the rarest quasars and the fundamental MBH-Mb„i g£ correlation in a hierarchical 
Universe. 

Subject headings: black hole physics — cosmology: theory, early Universe — galaxies: active, formation, evo- 
lution, starburst, high-redshift, ISM — methods: numerical — quasars: general, individual 
(SDSS Jl 148+5251) 



1. INTRODUCTION 

Quasars rank among the most luminous objects in the 
Universe and are believed to be p owered by SMBHs (e.g., 
ISalpeteri[l96l lLynden-Belllll969l) . They constrain the for- 
mation and evolution of galaxies and SMBHs throughout 
cosmic time. The s imilarity between t he cosmic star for- 
mation history (e.g., iMadau et ail 1 19961 iBunker et aHl2004t 
iBouwens et al.ll2004l) and the evolution of quasar abundances 
(e.g.. IShaver et al.l 119961) suggests an intriguing link be- 
tween galaxy formation and black hole growth. This is 
strengthened by tight correlations measured locally between 
the masses of the black holes and the global properties of 
the spheroid compo nents of their hosts, such as their lumi- 
nositi es and masses (Ma gorrian et al.lfl9 98; Mar coni & Hunt! 
2003), light concentration dGraham et al.ll2.001b, and velocity 
dispersions dFerrarese & Merrittll2000n Gebhardt et al.1 120001: 
iTrem aine et al l2002l) 

Distant, highly luminous quasars are important cosmologi- 
cal probes for studying the first generation of galaxies, the star 
formation history and metal enrichment in the early Universe, 
the growth of the first supermassive black holes (SMBHs), the 
role of feedback from quasars and black holes in galaxy evo- 
lution, and the epoch of reionization. The Sloan Digital Sky 
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Survey (SDSS. lYork et al.ll2000l) has contributed significantly 
to the discovery of high redshift quasars. Currently, there are 
over 1000 quasars know n at z > 4 and 1 2 at z > 6 (|Fan et al.l 
20011 12003112001 12006). As reviewed byEi(l2006), quasars 
at z ~ 6 are characterized by: (1) a low space density (~ 
10~ 9 Mpc~ 3 comoving); (2) high luminosities (absolute lu- 
minosity at rest-frame M1450 < -26), believed to be pow- 
ered by SMBHs of ~ lQ 9 Mp,\ ( 3) Gunn-Peterson absorption 
troughs dGunn & Petersonl 1965b in their spectra, which place 
these qua sars at the end of t he epo ch of reionization (e.g 



Fan et al 



Lidz et al. 



2001; Becker et al. 2001; Diorgovski et al. 2001 
2002; Songaila & Cowie 2002:" Wrriteet all 12003 



Sokasian et al.l l2003): and (4) a lack of evolution in the spec- 



tral e nergy distribution compared to lower-shift counterparts 
(e^. lElvis etai] [19941 iGlikman et al.l 120061 iRichards et ail 
2006), which demonstrates the existence of "mature" quasars 
at early times and comparable metal enrichment in quasars at 
all cosmic epochs. 

The most distant quasar known, SDSS Jl 148+5251 (here- 
after Jl 148+525 1), was discovered by SDSS at z = 6.42 
dFan et al.l l2003). It is extremely bright optically with M1450 
= - 27.8, and deep imaging surveys in both optical and ra- 
dio dCarilli et alJl2004t IWhite et"aLll200l IwTllott et alJl2005h 
show no sign of gravitational lensing or other companions 
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at the same redshift in the vicinity. Over the past few 
years, this quasar has been extensively st udied at many wave- 
leng ths. Near-infrared observations by IWillott et al.l d2003l) 
and iBarth et ail d2003l) imply a bolometric luminosity of 
^boi ~ 1O I4 L0 powered by accretio n onto a SMBH of m ass 1- 
5 x 10 9 Mf7y Radio observations by Bert oldi et all (l2003al) and 
ICarilli etalj d2004l) suggest that the host is a hyper-luminous 
far-infrared (FIR) galaxy, with Lfir ~ 10 13 L©, and these au- 
thors estimate a star formation rate of ~ 3 x 10 3 M©yr _1 by 
assuming most of the FIR luminosity comes from young stars. 
Emission from carbon monoxide (CO) has been detected 
dWalter et alJl2003t iBertoldi et al.ll2003bl: IWalter et al.ll2Q0l 
corresponding to a mass o f ~ 2 x 10 10 M < T). Dust has been seen 
by several groups (e.g.. IRobson et al.l 120041 IBertoldi et al.l 
l2003at ICarilli et alj |20o4 iBeelen et all l2006h with an esti- 
mated ma ss of ~ 5 x 10 8 MfD . In p artic ular, Spitzer obser- 
vations bv lCharmandaris et all (120041) and lHines et all d2006l) 
indicate that the dust is heated by an active gala ctic nucleus 
(AGN ). Furthermore, the detection of iron by Bart h et al 
( 20031) . the carbon fine structure lin e rCIIl bv iMaiohno et al 



( 2005b and excess OI absorption by [Becker et al.l d2006b indi 
cate near-solar metallicity in the quasar host. 

These observations raise many fundamental questions for 
models of quasar and galaxy formation: Where do such high- 
redshift, luminous quasars originate? How do they form? 
What are the mechanisms and physical conditions for SMBH 
growth? And, do these quasar hosts obey the same SMBH- 
host correlations as observed in the local Universe? 

Interpretations of various observations of J 1148+5251 
have painted a rather conflicting picture for the forma- 
tion site of the quasar halo and the SMBH-host relation- 
ship. The low abundance of these quasars leads to the 
view that they are hosted by massive halos (> 10 I3 M©) 
in the rarest de nsity peaks of the dark matter distribution 
dFan et al.l 120031) . However, it has been argued, based on 
the lack of companion galaxies in the field, that this quasar 
may reside i n a much lower mass halo in a le ss over- 
dense re£ionJCarilhetap|200l iMUottitiDIIOOll). More- 
over, Walt er et all d2004l) suggest, based on the dynamical 
mass estimate from CO measurements, that J 1 148+5251 con- 
tains a small stellar spheroid, and that the SMBH may have 
largely formed before the host galaxy. However, the de- 
tection of metal line s dWalter et al.l l2004t IBarth et al l l200l 
Maioli no et"aT1l2005l). along with dust (IBertoldi et alJl2003at 
Carilli et all 120041; IRobson et all 12004 ICharmandaris et all 
2004tlHines et al.ll2006l;lBeelen et al.ll2006l) . indicates that the 
interstellar medium (ISM) of J 1148+5251 was significantly 
enriched by heavy elements produced through massive star 
formation at rates of ~ 10 3 MQyr _1 occurring as early as 
z > 10, and that large stellar bulges form before accreting 
SMBHs undergo luminous quasar activity. 

In an expanding Universe that is dominated by cold dark 
matter and is accelerated by dark energy, the A-cold dark 
matter (ACDM) cosmology, the leading theoretical model for 
structure formation, assumes that structure grows from weak 
density fluctuations amplified by gravity, with small objects 
collapsing first and subsequently merging to form progres- 
sively more massive ones, a process know n as "hierarchical 
assembly" (e.g., see lBarkana & Loebll2001l for a review). The 
formation of galaxies and quasars is therefore determined by 
the abundance of dark matter halos; i.e., the number den- 
sity of halos as a function of mass and redshift. The most 
widely used analytic model for the halo ma ss function was 
first developed by iPress & Schechted d 19741) (hereafter PS), 



which is based on Gaussian density perturbations, linear grav- 
itational growth, and spherica l collapse of dark matte r. Us- 
ing the PS mass functions, Efstathiou & Rees (1988) stud- 
ied the abundance of rare objects, such as luminous quasars 
at high redshifts. These authors predicted a sharp "cut-off" 
of the quasar population at z ~ 5. However, while the ini- 
tial, linear growth of density perturbations can be calculated 
analytically, the gravitational collapse and subsequent hierar- 
chical build-up of structure is a highly nonlinear process that 
can be followed only through numerical simulat ions. It has 
been s hown by previous numerical studies (e.gJJenkins et al.l 
120011: IShem & Tormenl 120021: ISpringel & Hernquistl l2003b). 
and more recently by the state- of-the-art Millennium Simu- 
lation by ISpringel et al.l d2005d) . that the PS formula under- 
estimates the abundance of high-mass halos by up to an order 
of magnitude. Therefore, whether or not rare quasars such as 
Jl 148+5251 can form in the ACDM cosmology remains an 
open question and an important test of the theory. 

To date, only a limited number of analytical or semi- 
analytic models have addressed the early formation of a 
10 9 M P) SMBH at z ~ 6 dHaiman & Loebl200UlHaimanl2004t 
lYoo & Miralda-EscudeH2004tlVolonteri & Reesll2005l) . These 
approaches use merger trees of dark matter halos generated 
using the PS theory, and assume a black hole accretion rate at 
or above the Eddington limit. However, as discussed above, 
the PS-based approach may be inaccurate. Moreover, it is not 
clear whether sufficient physical conditions for such large ac- 
cretion rates exist in quasar systems as the hydrodynamics of 
the large-scale gas flow and feedback from black holes have 
not been incorporated in earlier modeling. 

It is believed that the growth of SMBHs is c losely 
linked to galaxy formation (e.g.. [M agorrian e t al.l 119981; 



Ferrare se & Merritt 2000; Geb hardt et al . 2000; Gr aham et all 
20011 iTYemaine et al.l 120021; iMarconi & Hunt! 120031 
Haimanl l2004t iKazantzidis et al.l 120051: iLietaiT [2006, 



and that th e grow t h is self-regula t ed by feedback (e.g. , 
Silk & Reesl [19981; iHaehnelt et al ' 



1999 
2005 



. 19981; iF abian 

Kind 120031; IWvithe & Loebl 120031; [Pi Matteo et al 
Spring eTet al.l l2005bl; ISazonov et al.1 120051; iMurrav et al.1 
120051: IWvithe & Loebl 12005b . Remarkably, self-regulated 
models with SMBH feedback in the form of thermal energy 
coupled to the ambient gas have been demonstrated to suc- 
cessfully reproduce many observations of galaxies, including 
the Mb h-c relation dPi Matteo et all [20051 IRobertson et al 
2006b), galaxy colors (Sprin gel et al.1 2005al; IHopkins et al. 
2006d) . X-r ay gas emission dCox et al.l l2006bb . elliptical 
kinematics dCox et al.l 120 06c) and the fund amental plane 
dRobertson et al.l l2006ab . quasar p roperties dHopkins et al.l 
luminosity f unctions (H opkins etaL[_ I2005bld. 



2006b), and populations (Hopkins et al. 2006a e; J), and the 
lumin osity function of low-level AGN (Hopkins & Hernquist 
2006). Furthermore, these simulations of binary mergers 
identify a plausible, merger-driven formation mecha- 
nism for massive black holes and luminous quasars (e.g., 
| Pi Matteo etail 120051: IHopkins et al.ll2006ab IRobertson et all 
20063). 

Here, we present a model that accounts for the SMBH 
growth, quasar activity and host galaxy properties of the most 
distant quasar known. In our scenario, the quasar and its host 
galaxy form in a massive halo that originates from a rare den- 
sity peak in the standard ACPM paradigm, and they grow 
hierarchically through multiple gas-rich mergers, supporting 
an average star formation rate of ~ 10 3 M Q yr" 1 that peaks at 
z ~ 8.5. Once the progenitors undergo sufficient dynamical 
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friction to coalesce, the multiple SMBHs from the progeni- 
tor galaxies merge and exponentially increase their mass and 
feedback energy via accretion. At z ~ 6.5 when the SMBH 
mass exceeds 10 9 M Q , black hole accretion drives a suffi- 
ciently energetic wind to clear obscuring material from the 
central regions of the system and powers an optically lumi- 
nous quasar similar to Jl 148+5251. We have devised a set 
of novel multi-scale simulations, which include cosmological 
N-body calculations on large scales and hydrodynamic simu- 
lations of galaxy mergers on galactic scales, coupled with the 
self-regulated growth of SMBHs, enabling us to follow galaxy 
assembly and quasar formation at z ~ 6. 

This paper is organized as follows. In § 2, we describe 
our computational methods and models, which includes a set 
of large scale cosmological N-body simulations, and hydro- 
dynamical galaxy mergers along the merging history of the 
quasar halo. In § 3, we present the formation and evolution 
of the quasar and its host galaxy, including the assembly of 
the galaxy progenitors, star formation, and SMBH growth, 
as well as the SMBH-host correlations, and properties of the 
quasar such as luminosities and lifetimes. We discuss feed- 
back from starburst-driven winds, quasar abundances for cos- 
mological models with different parameters, the implication 
of black hole mergers, and galaxies in the epoch of reioniza- 
tion in § 4, and summarize in § 5. 

2. METHODOLOGY 

Rare, high-redshift quasars originate in highly overdense 
regions in the initial dark matter density distribution and grow 
through hierarchical mergers, as predicted by the ACDM 
theory. Simulations of high-redshift quasar formation must 
model a large cosmological volume to accommodate the low 
abundance of this population, have a large dynamic range to 
follow the hierarchical build-up of the quasar hosts, and in- 
clude the hydrodynamics of the gas flows in g alaxy mergers. 
The cu tting-edge Millennium Simulation by ISpringel et alJ 
(2005c) follows structure formation in a box with side 
length of 500 /T 1 Mpc using 2160 3 dark matter particles. It 
reproduces the large -scale galaxy distribution as observed 
( Sprin gel et al.l l2006). and identifies an early quasar halo can- 
didate at z = 6.2 which ends up in the richest cluster at the 
present day. However, even such a large dynamic range still 
falls short of being able to follow the formation and evolution 
of the rarest quasars observed at the highest redshifts. More- 
over, in order to address the properties of quasars and host 
galaxies, gasdynamics and physical processes related to star 
formation and black hole growth must be included. To satisfy 
these requirements, we have performed a set of novel multi- 
scale simulations that enable us to resolve individual mergers 
on galactic scales and retain the context of large-scale struc- 
ture formation, as well as the evolution of black holes and 
stars. 

First, we perform a coarse dark matter simulation in a vol- 
ume of 1 /T 3 Gpc 3 designed to accommodate the low number 
density of z ~ 6 quasars. The largest halo at z = 0, within 
which the desce ndants of early, lumi nous quasars are as- 
sumed to reside (Springe l et al.l 1200 5c). is then selected for 
resimulation with higher re solution us i ng a m ulti-grid zoom- 
in technique developed by iGao eT al. (2005). The merging 
history of the largest halo at z ~ 6, which has reached a 
mass of ~ 5.4 x 10 12 /r' M© through 7 major (mass ratio < 
5:1) mergers between redshifts 14.4 and 6.5, is extracted. 
These major mergers are again re-simulated hydrodynami- 
cally using galaxy models scaled appropriately for redshift 



dRobertson et al. 2006b) and adjusted to account for mass ac- 
cretion through minor mergers . The simulations include pre - 
scriptions for star formation ( Springel & Hernquist 2003a), 
and for SMBH grow th and feedback dDi Matteo et alj|2005t 
Spri ngel et al.ll2005bl) . as described below. 

2. 1 . Code and Parameters 

Our multi-scale simulations were performed using the par- 
allel, N-body/Smoothed P article Hydrod ynamics (SPH) code 
GADGET2 developed by ISpringel! (120051) that is well tested 
in a wide range of applications from large scale structure for- 
mation to star formation. For the computation of gravita tional 
forces, the code uses the "Tr eePM" method (|Xu[l 19951) that 
combines a "tree" algorithm dBarnes & Hutlfl986l) for short- 
range forces and a Fourier t ransform particle-mesh method 
dHockney & Eastwo odll 19811) for Ion g-range forces. The PM- 
method works efficiently in large scale cosmological simula- 
tions, while the tree-method provides accurate forces for the 
large dynamic range of galaxy merger simulations. 

GADGET 2 implements the entropy-c onserving formula- 
tion of SPH (Spring el & Hernquistll2002l) with a daptive par- 
ticle smoothing, as in Hernquist & Katz (1989). Radiative 
cooling and heating processes are calculated as s uming col- 
lisional ionization equilibrium dKatz et al.l 119961: iDave et all 
1999). Star formation is modeled in a multi-ph ase ISM, 
with a rate that follow s the Schmidt-Kennicutt Law (ISchmidtl 
119591: lKennicuttlfT998T) . Feedback from supernovae is cap- 
tured through a multi-phase model of th e ISM by an effective 
equatio n of state for star-forming gas ( Springel & Hernquist 
l2003al) . A prescription for supermassive black hole growth 
and feedback is also included, where black holes are rep- 
resented by collisionless "sink" particles that interact gravi- 
tationally with other components and accrete gas from their 
surroundings. The accretion rate is estimated from the local 
gas density and sound speed using a spherical Bondi-Hoyle 
( lBondi&Hovle1ll944t iBondill 19521) model that is limited by 
the Eddington rate. Feedback from black hole accretion is 
modeled as the rmal energy in j ected i nto s urrounding gas, 
as des cribed in Springel et al.l d200 5b) and Di Matte o et alJ 
J2005I 

The simulations presented in this paper adopt the ACDM 
model with cosmological parameters chosen according to 
the first y ear Wilkinson Micr owave Anisotropy Probe data 
(WMAPl. lSpergel et alJl2003h . (O m , O b , fi A) h, n s , cr 8 )= (0.3, 
0.04, 0.7, 0.7, 1, 0.9). Here, Q m is the total matter density 
in units of the critical density for closure, p cr ; t = 3>Hq/{%ttG). 
Similarly, 51b and ft\ denote the densities of baryons and 
dark energy at the present day. The Hubble constant is pa- 
rameterized as Ho = 100/ikms" 1 Mpc -1 , while er 8 is the root- 
mean-squared linear mass fluctuation within a sphere of ra- 
dius 8 h~ l Mpc extrapolated to z = 0. We have also done the 
same cosmolog ical simulations wi th WMAP third year re- 
sults (WMAP3,Hisrg§Let al. 2006), (Q m , ft h , ft A , h, n s , cr 8 )= 
(0.236, 0.042, 0.759, 0.732, 0.95, 0.74) for comparison. 

2.2. Cosmological Simulations 

The quasars at z ~ 6 have an extremely low comoving space 
density, n ~ 10~ 9 Mpc~ 3 , and are believed t o reside in mas - 
sive dark matter halos with M > 10 13 M Q dFan et al.l l2003h . 
Cosmological simulations of quasar formation must therefore 
model a volume of ~ 1 h~ 3 Gpc 3 to account for the rarity of 
such objects. However, in order to resolve a 1O 13 M halo at 
z ~ 6 in a cosmological simulation with uniform resolution, 
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FIG. 1 . — Snapshots from a cosmological simulation run with WMAP1 parameters. The images show projected density of dark matter in x-y (left column) and 
x-z (right column) planes, the red dot represents the center-of-mass of the quasar halo, which is the largest halo at both z = and z = 6. The top panels show 
the coarse run at z = 0. The middle and bottom panels show the zoom-in run at z = 6.06 and z = 0, respectively, the number at the lower-left corner indicates the 
number of groups identified at that redshift. 
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a dark matter particle mass at least as small as 10 A Mq 
and particle numbers of > 10 9 are required. Tracking the 
merger history of such halos requires ~ 2 orders of magnitude 
higher resolution and would be computationally prohibitive 
with standard techniques. 

We achieve the mass resolution requirements for the merger 
history of a 10 13 M© halo at z ~ 6 by means of a two-step 
re-simulation. First, coarse dark matter cosmological sim- 
ulations are performed to identify a candidate halo for the 
quasar host. A cubic volume L\, ox = 1 h~ l Gpc on a side is 
simulated with 400 3 particles, achieving mass and force res- 
olutions of md m ~ 1.3 x 10 12 h~ l Mq and e ~ 125/r'kpc (co- 
moving), respectively. To generate the initial conditions, we 
use th e Boltzmann code CMBFAST by Seliak & Zaldarriaga 
( 1996) to compute a linear theory power spectrum for our cho- 
sen cosmology. A random realization of the power spectrum 
is constructed in Fourier space, sampling modes in a sphere 
up to the Nyquist frequency of the mesh. The particle dis- 
tribution is evolved forward in time to z = from its initial 
displacement at z = 30 determined using the Zel'dovich ap- 
proximation. 

At the end of the simulation, halos are identified us- 
ing a "friends-of -friends" (FOF) group-finding algorithm 
dDavis et al.l 19851) with a fixed comoving linking length equal 
to 0.2 times the mean dark matter interp article separation and 

mini mum of 32 particles per group ( Sprin gel & Hern quist 



2003b). The mean overdensity of the groups corresponds 
approximately to the expected density of virialized halos 
(Sprin gel et al.ll2005ch . From the more than 126000 groups 
identified in the 1 h~ 3 Gpc 3 volume at z = the largest halo 
with M(z = 0) ~ 3.6 x 1O 15 /i _1 M is selected as a candidate 
halo for modeling the formation of a qu asar at z = 6.5. 

A multi-g r id tech nique developed by iGao et all d2005l) and 
iPower etaO (120031) is used to " zoom in" with high resolu- 
tion on the selected halo region which has an effective side 
length of Lbox ~ 50/T 1 Mpc. Large-scale tidal forces are cap- 
tured by binning exterior particles into cells according to their 
distance from the high-resolution region. To ensure proper 
treatment of small-scale structure, the initial displacements of 
the high resolution particles are calculated assuming a higher 
initial redshift of z = 69 and normalized to a% at z = 0. The 
re-simulation uses w 350 3 particles, with w 340 3 particles in- 
side the high-resolution region. With this technique, the mass 
resolution increases by almost four orders of magnitude to 
nidm ~ 2.8 x 10 8 /z -1 Mq while the spatial resolution reaches 
e ~ 5 h~ l kpc. 

Figure Q] shows snapshots of both the coarse and high- 
resolution zoom-in runs that locate the quasar halo candi- 
date. In the coarse run, the "cosmic web" is clearly seen, al- 
though the distribution appears nearly homogeneous on such 
large scales. In the zoom-in run, filamentary structures are 
prominent. Dark matter collapses along the filaments, and the 
largest halo forms in the deepest potential wells at the inter- 
sections of the filaments. The high resolution of the zoom-in 
run enables the identification of more halos with lower masses 
both at z = and at high redshifts as early as z ~ 17, which 
is sufficient to identify the halo progenitors of the candidate 
quasar at z ~ 6. It appears that the halo progenitor of the 
largest one at the present day is also the most massive halo at 
z ~ 6, when it reaches a mass ofM w 5.4 x 1O 12 /! _1 M , mak- 
ing it a plausible candidate for hosting a rare z ~ 6.4 quasar. 

2.3. Halo Mass Functions with Different Cosmological 
Parameters 



The impact that variations in the cosmological parame- 
ters can have on large-scale structure formation can be un- 
derstood from the theoretica l mass function of halos, as de- 
rived by iPress & Schechterl (1 19741) and later developed by 
lLacey & Cold (|1993|). The comoving number density dn of 
halos of mass between M and M+dM can be described as, 



dn 
dM 



2 pa_ S c (z) 
7T M 2 cr(M) 



d lner 



dlnM 



exp 



S c (z) 2 



2a 2 (M) 



(1) 



where po is the local mean mass density, 8 c (z) is the critical 
density of collapse at redshift z linearly extrapolated to the 
present day, and cr(M) is the mass variance, which is a func- 
tion of the power spectrum P(k) with wavenumber k and the 
window function w(k), a 2 (M) = J °° P(k)w 2 (k)d 2 k. The 
abundance of halos depends on the two functions a(M) and 
#crit(z), each of which involves the cosmological parameters, 
in particular a%, Vt m and f^A. These parameters determine the 
formation epoch of a halo and its mass. 

The recen t ly-rele ased third year WMAP3 results 
( Sperg eTet alj 120061) have lower values of a%, n s and 
Q m , compared to WMAP1 dSpergel et al.ll2003l) . The smaller 
(is from WMAP3 would lower the amplitude of the power 
spectrum, which in turn reduces er(M). Furthermore, a 
smaller H m would reduce S c (z) and hence delay halo for- 
mation. So, compared to WMAP1, at a given redshift the 
WMAP3 parameters would yield a lower abundance of halos 
with mass M^\ >M^,, where M* is the halo mass corre- 
sponding to the characteristic l uminosity in the Schechter 
luminosity function for galaxies dSchechterll 19761) ; while for 
Mhaio < M*, it predicts a slightly larger halo abundance. 

To test the sensitivity of our model to the new WMAP re- 
sults, we have performed the same set of cosmological sim- 
ulations with param eters from the WMAP3 measurements 
(Sperg eTet al.l20 06). We find that indeed the changes implied 
by the new parameters primarily affect the formation time and 
the mass of the candidate quasar halo. For the same random 
phases in the initial conditions, the location of the most mas- 
sive halo at z = remains the same in both the WMAP1 and 
WMAP3 runs, except that its mass is reduced by a factor of 
^1.6 for the WMAP3 parameters. Similarly, the mass of 
the largest halo at z ~ 6 is altered by roughly the same fac- 
tor. Other notable changes include: (1) the formation epoch 
of the first halo is shifted from z ~ 16.8 in the WMAP1 run 
to z ~ 14.4 in the WMAP3 run, and (2) the merging history 
of the largest halo at z ~ 6 moves to lower redshifts in the 
WMAP3 run, but the number of major mergers remains the 
same. 

Figure [2] shows the halo mass functions from differ- 
ent cosmological simula tions. The PS mass function 
dPress & Schechterl 1 1974b . as well as the one corrected 
to match numerical simulations by ISheth & Tormenl d2002l) 
(ST), are also shown for comparison. One important feature in 
this figure is that the coarse runs agree well with the ST mass 
function, but show a larger comoving density at the high mass 
end than that predicted by the PS theory. Our results show 
that the PS formula underestimates the abundance of high- 
mass halos by nearly an order of magnitude at z = 0, and the 
discrepancy between the PS calculation and numerical simu- 
lations becomes l arger at higher redshifts, confirming previ- 
ous findings (e.gJJenkins et alJl200lUSheth & Tormenll2002l: 
ISpringel & Hernquist! l2003bl: ISpringel et all l2005d) . This 
may explain why previous models using the PS formula to 
study the abundance of luminous quasars, which presumably 
form in massive halos, under-predicted the number of bright 





FIG. 2. — Halo mass functions from cosmological simulations with parameters from WMAP1 (left) and WMAP3 (right), respectively, and with different levels 
of resolution corresponding to our coarse (top panels) and zoom-in (bottom panels) runs. The colored symbols indicate different redshifts, while the error bar 
shows the Poisson error VN. The mass functions from Press & Schechter 1 1974) (PS, dotted line) and Sheth & Tormen (2002) (ST, solid line) are also shown for 
comparison. Please note that in the bottom panels, the analytical curves (ST and PS) apply only to a random region, they are not suitable for a highly overdense 
region where the most massive halos reside in the zoom-in box. So the high-mass end of the simulated mass function deviates significantly from the prediction, 
see text for more discussion. 



quasars at z > 5 (e.g., Efstathiou & Rees 1988). Furthermore, 
these results also suggest that the commonly used analytical 
merger tree generated using the PS formula may not be suit- 
able to study quasar formation at high redshifts. 

There are two clear "shifts" of the mass function caused 
by resolution and cosmological parameters. Those from runs 
with higher resolution extend to higher redshifts, and at the 
same redshift, the WMAP1 runs produce more massive halos 
than the WMAP3 ones. As shown in Figure [2] the coarse runs 
produce mass functions only up to z ~ 3 owing to limited mass 
resolution, while the zoom-in runs can produce quite reason- 
able mass functions as early as z ~ 14. Because the zoom-in 
runs were deliberately centered on the highest density peak 
of the 1 h~ 3 Gpc 3 box, they each contain a very massive halo 
(M > 10 15 M Q atz = 0)by construction. This explains why the 
highest mass bin (which contains only one halo in this case) 
is ~ 2 orders of magnitude larger than the theoretical curves 
(ST, PS), which apply only to a random region that has a much 



lower density fluctuation. 

To summarize, at a given redshift, runs with the WMAP3 
parameters yield slightly less massive halos than ones per- 
formed with the WMAP1 values. Or, to put it differently, 
objects in the WMAP3 cosmology will have masses similar 
to those for WMAP1, but at slightly later times (i.e. lower 
redshifts). In what follows, we are primarily concerned with 
investigating the plausibility of forming z ~ 6 quasars through 
the self-regulated growth of SMBHs in hierarchical mergers, 
rather than precisely reproducing the properties of an individ- 
ual quasar at a given redshift, such as J 1148+5251. Most of 
our results are therefore based on runs with the WMAP1 pa- 
rameters, to ease comparison with previous numerical work. 
If it were firmly established that e.g. as is in reality smaller 
than its WMAP1 value, then a more exact match to a particu- 
lar quasar could presumably be obtained by considering larger 
simulation volumes and identifying a suitable candidate host 
that is slightly rarer than the one we have chosen to focus on 
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FIG. 3. — Schematic merging history of the largest halo at z = traced by 
mergers at different redshifts with mass ratio < 10 : 1, which is defined as 
the mass ratio between the halo and progenitor at a given time. Each of the 
progenitors joins in this big merging event at a given redshift, interacts with 
the system, and subsequently merges with others at later times. The quasar 
host at z ~ 6 is built up by seven successive major mergers of progenitors Gl, 
G2, ... G8 from z ~ 14.4 to z. ~ 6.0, as illustrated by the red lines in this plot. 
The first interaction between Gl and G2 takes place at z ~ 14.4, then G3 and 
G4 join in the system at z, ~ 13, followed by G5, G6, G7 and G8 at later times 
(see text for more details). The timeline of these events, the mass and other 
properties of these progenitors are listed in Table 1 . 



here. 



2.4. Merger Tree Construction 



To follow the hierarchical mass assembly of the host galaxy 
over cosmic time, the merger tree of the halo is extracted from 
the cosmological simulation. This tree provides key infor- 
mation for computing the physical properties of the progen- 
itor galaxy population. While the merger history of the halo 
includes a spectrum of progenitor masses, the most massive 
progenitors contribute the majority of the halo mass over the 
redshift range considered. We trace the merger history of the 
most massive progenitor at each redshift by using particle tags 
to identify progenitor systems at earlier redshifts in the sim- 
ulation. Groups that contribute at least 10% of the halo mass 
at a given time step are considered as the progenitors of the 
halo and are recorded. The procedure is repeated until the last 
progenitor is reached, producing the merging history. 

Figure [3] illustrates the merging history of the largest halo 
at z = in our cosmological simulation, which has a mass 
of ~ 3.6 x 10 15 /z -1 M Q . It is also the largest one at z ~ 6 
with a mass of ~ 5.4 x 10 12 h M©. This schematic plot out- 
lines the redshift of merger event, and the mass ratio of the 



halo to its galaxy progenitors at a given time. It shows that 
this halo grows rapidly through hierarchical mergers early on, 
with seven major mergers (mass ratio of the merging pairs 
< 5 : 1) from z ~ 14.4 to z ~ 8.5 that build up a substantial 
fraction of the halo mass at z ~ 6. 

In modeling the development of a z ~ 6 quasar, we are pri- 
marily interested in "major" mergers, where the mass ratio 
of the merging galaxies is not too far from unity, for sev- 
eral reasons. First, it is believed that major mergers play the 
most important role in the formation and evolution of mas- 
alaxies (e.g.. ISanders & Mirabellll996t IScoville et al 



sive 



2000UVeilleux et al. 2002; Conselice et al. 2003; Dasy ra et al 



2006). Second, and of greater concern to us in this paper, in 
our picture for quasar fueling, gas in a rotationally supported 
disk loses angular momentum through gravitational torques 
excited by tidal forces in a merger, driving the growth of su- 
permassive black holes. This process operates most effec- 
tively in a major merger because the tidal deformation of each 
galaxy is si g nifica nt in such an event (Ba rnes & Her nquist 
1 1 99 lL 1 19921 119961) . Collisions involving galaxies with a 
mass ratio as large as 10 : 1 can induce ga s inflows in disks 
dHernquistlll989HHernquist & Mihosll 19951) . but only for lim- 
ited orbital configurations. For these reasons, we focus on 
mergers from the merger tree having a mass ratio < 5 : 1, as 
outlined by the red color in Figure [3] 

In the resimulation of the merger tree as described in § [ 
we take into account mass accretion of the halo by adding 
mass proportionally to each of the eight progenitor galaxies 
in the major mergers. This approach preserves the progenitor 
mass ratios and approximately preserves the dynamics of the 
major mergers (Dubinski 1998|). 

2.5. Simulations of Galaxy Mergers Along the Tree 

In order to model the formation and evolution, and prop- 
erties of the quasar candidate, the merger tree is then re- 
simulated hydrodynamically with galaxy models that consist 
of an extended dark matter halo, a rotationally supported, ex- 
ponential disk of gas and stars and a central supermassive 
black hole. We follow the evolution of the system built up 
by seven major mergers hierarchically from z ~ 14.4 to z ~ 6, 
as shown in Figure [3] Technically, this is a series of succes- 
sive merger simulations. The first simulation includes Gl and 
G2 interacting at z ~ 14.4. It stops at z ~ 13 and a new galaxy 
G3 is added into the system. During this process, all the dy- 
namical properties of the pre-existing system (e.g. Gl and 
G2 in this case) are preserved, while G3 is added based on 
its properties and orbital parameters derived from cosmolog- 
ical simulations. Then a second merger simulation with Gl, 
G2 and G3 starts. Such a procedure is repeated until all the 
progenitors enter the system. In the end, the simulation in- 
cludes all eight galaxies. Eventually all these galaxies and 
black holes merge together. The duration of each merger sim- 
ulation is determined by the merger tree. The redshift at which 
each progenitor galaxy enters the merger tree, the properties 
of each progenitor galaxy, and the numerical parameters of the 
merger simulations are listed in Table 1 . Below we describe 
the specification of these parameters. 

2.5.1 . Galaxy Models 

The structure of the galaxy models is motivated from lead- 
ing theories of dissipational d isk galaxy formation in CDM 
cosmologies that, as shown by lMo et aD (11998b . are success- 
ful in reproducing the observed properties of both present-day 
disk galaxies and damped Lya absorbers in quasar spectra at 
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TABLE 1 

Progenitor Properties and Numerical Parameters 



Gnl<ixy a 


Z 


M ■ c 


v ■ d 


f e 
Jgas 




K v 


ff„ h 






[10 10 /i-' M Q ] 


[kms- 1 ] 




[io 5 r'M ] 


[fr 1 kpc] 


[hr 1 kpc] 


Gl 


14.4 


6.3 


234.1 


1.0 


0.15 






G2 


14.4 


5.3 


220.3 


1.0 


0.15 


0.2 


7.1 


G3 


13.0 


15.0 


297.8 


1.0 


0.51 


0.2 


8.5 


G4 


13.0 


17.7 


314.6 


1.0 


0.51 


0.3 


10.7 


G5 


10.5 


49.1 


401.0 


1.0 


6.56 


0.4 


11.3 


G6 


9.4 


79.6 


448.6 


0.9 


23.3 


0.5 


18.2 


G7 


8.5 


160.0 


540.4 


0.9 


89.2 


0.7 


25.2 


G8 


8.5 


207.7 


589.5 


0.9 


89.2 


1.0 


34.5 



a Name of galaxy progenitor. Gl is the halo at z = 14. 4. b Redshift at which the progenitor enters 
the merger tree. c Virial mass, assuming overdensity A = 200. d Virial velocity, assuming overden- 
sity A = 200. c Gas fraction of the progenitor baryonic mass/ Progenitor black hole mass at the 
merger redshift. 8 Pericentric distance of the incoming progenitor to the center-of-mass of the exist- 
ing system. h Initial distance of the incoming progenitor to the existing system. 



high redshift. The initial galaxy models are con structed in 
dynamical equilibrium using a well-tested method dHer nquist 
19931; [Springel & Whitdll99l ISpringell2000t ISpringel et all 
2005b). A halo is identified with a virial mass M wu and a virial 
radius /? v j r within which the overdensity A = po/Pcrit = 200, 
where po and p a i t are the mean and critical density, respec- 
tively. The density profile of the da rk matter halo follows 
a Hernquist profile (Hernquist! 119901). scaled to match that 
found in co smological si mulations~i lNavarro eFaIlfT997h . as 
described in Springel et al.ld2005bl) : 



M v 



(2) 



27r r(r+a) 3 ' 

where a is a parameter that relates the Hernquist ( 1990) profile 
parameters to the appropriate NFW halo scale-length R s and 
concentration C v i r (C v i r = R W u/R s ), 



a = R s v/2[ln( 1 + C vll -) - C vll -/( 1 + C vil -)] . 



(3) 



T he exponen t ial dis k of stars and gas are then constructed 
as in iHernquistl (119931) and lSpringel et all (l2005bl) . The prop- 
erties of the galaxy, including the virial mass M v i r , virial ra- 
dius R vlI and halo concentrati on C wlr are scaled appro priately 
with redshift, as described in Robertson et al. (2006b). In par- 
ticular, for a progenitor with virial v elocity V v i r at red shift z, 
M vir and R VK are cal culated followinglMo et all (fl998). while 
Cyii is adopted from lBullock et all (120011) . as briefly outlined 
below: 



l0GH(z) ' 

Vvir 



(4) 

l0H(z) ' (5) 
H(z) =H [n A + (1 - Q A - O m )(l +z) 2 + n m ( 1 +z) 3 ] 1/2 , (6) 

-0.13 

(1+zT 1 . 



M v 



Mn 



(7) 



where G is the gravitational constant, and Mo ~ 8 x 
10 12 h Mq is the linear collapse mass at the present epoch. 

We assume a baryon fraction of ft, = 0.15 f or these high- 
redsh ift galaxies based on the WMAP1 result ( Sperg el et al.l 
2003). The gas fraction of each progenitor is extrapolated 
from the results of sem i-analytic models of galaxy formation 
dSomerville et alj|2001l) . with 100% gas disks at z > 10 and 



90% at 10 > z > 8. The multiphase ISM is envisioned to con- 
sist of cold clouds embedded in a hot, tenuous gas in pressure 
equilibriu m. Stars form ou t of the cold clouds by gravitational 
instability dLi et al.l 12005) wit h a rate that is proportional to 
the surface den sity of the gas (ISchmidtl 1 959t iKennicuttl 1 9981: 
lLietal.ll2006bl) . 

In the adopted ISM model for the gas, the equation of state 
(EOS) is controlled by a parameter g E os that linearly interpo- 
lates between isothermal gas (qeos = 0) and a strongly pres- 
surized multiphase ISM (geos = 1)- This EOS describes the 
dynamics of star-forming gas and accounts for the conse- 
quences of stellar feedback on galactic scales, and enables 
us to const ruct equilibrium disk models even with large ga s 
fractions dRobertson et al. 2004; Spring el & Hernq uist 2005). 
Supernova feedback is modeled through thermal energy in- 
put into surrounding gas and can help evaporate the cold 
clouds to replenish the hot phase. For the simulation pre- 
sented here a value of geos = 0.5 is used, but test simulations 
using <7eos = 0.25- 1.0 produce average star formation and 
black hole accretion rates that converge to within 15%. 

2.5.2. Black Hole Accretion and Feedback 

The supermassive black holes are represented by collision- 
less "sink" particles. They interact with other particles gravi- 
tationally, and accrete the gas. Accretion of gas onto the black 
holes is modeled us ing a Bondi-Hoyle-Lyttleton parameteri- 
zation dBondilll952l;lBondi & Hovldll944tlHovle & Lvttletonl 
1 19411) . in which the black holes accrete spherically from 
a stationary, uniform distrib ution of gas, as d escribed in 
Di Matte o et all d2005l) and lSpringel et alj d2005bl) : 



Mr = 



4naG 2 Ml H p 



(8) 



(c 2 + v 2)3/2 ' 

where Mrh is the black hole mass, p and c s are the density 
and sound speed of the gas, respectively, a is a dimensionless 
parameter of order unity, and v is the velocity of the black hole 
relative to the gas. 

We assume the accretion has an upper limit by the Edding- 
ton rate, 

47rGM BH OT p 



M Edd 



(9) 



e r CT T c 

where m p is the proton mass, aj is the Thomson cross- 
section, and e r is the radiative efficiency. The latter deter- 
mines the conversion efficiency of mass accretion into en- 
ergy released as radiated luminosity. We adopt a fixed value 
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of e r = 0.1, which is the m ean value for radiatively efficient 
Shaku ra~& Sunyaevl (1 19731) accretion onto a Schwarzschild 
black hole. In the simulations, the accretion rate is then the 
minimum of these two rates, Mbh = min(MEdd , Mb)- 

The feedback from the black holes is associated with the 
mass accretion. We assume that a small fraction (~ 5%) of 
the radiated energy couples to the surrounding gas isotropi- 
cally as feedback in form of thermal energy. This fraction 
is a free param eter, determined by m atching the observed 
Mbh-c relation (Di Matteo et al.l2005l). For m ore discussions 
on this prescription, see Hopki ns et al.ld2006al) . This feedback 
scheme self-regulates the growth of the black hole, and has 
been demonstrated to successfully reproduce many observed 
properties of elliptical galaxies, as mentioned earlier. 

2.5.3. Black Hole Seeds 

To grow a black hole up to 10 9 M Q in less than 
800 million years, a wide range in seed mas ses, from 
10M W to 10 6 M(7), have b een suggested (e.g.. ICarr et alJ 
1984 lLoeb & Rasiol 11994 iBromm & Loebl l2003b iHaimanl 
2004 lYoo & Mi ralda-Escude 2004; Volonteri & Reesj |2005t 
Begelman et al. 2006). The formation of the black hole seeds 
remains an open questio n, and seve r al sce narios have been 
proposed. In particular, Fry er et all (1200 ll) show that rapid 
collapse of massive PopIII stars due to pair instability c ould 
produce black hole of ~ 10 2 M Q ; Bromm & Loeb (2003]) sug- 
gest that hot and dense gas clump may collapse monolithicly 
to form a massive black hole of ~ 1O 6 M0 in metal-free galax- 
ies wit h a virial temperature of 10 4 K; while Be gelman et al.l 
(2006) propose that ~ 2OM black holes could form by di- 
rect collapse of self-gravitating gas due to global instabili- 
ties in protogalactic halos, they then grow to 10 4 ~ 6 M Q with 
super-Eddington accretion. We adopt the pic ture where black 
hole seeds are the remnants of the first stars (|Abel et al.ll2002t 
Bromm & Larson! 12004 iTan & McKeel 12004 lYoshida et al. 
2006t iGao et all 120061) . The remnant black hole mass is 
currently uncertain and widely debated. Recent theory of 
PopIII star formation predicts a mass range of ~ 30-500 M Q , 
but there are two regimes whe re a SMBH could form, ei- 
ther < lOOAfo or > 260 M (iHeger et all 120031, see also 
Yoshida et al. 2006 for a recent discussion). We have tested 
the seed mass in the range of 1 00-300 M Q and find that the ex- 
ponential growth of the black holes during the merger makes 
our results insensitive to the choice in that range. We there- 
fore assume that the black hole seed starts with an initial mass 
of 200 M Q after the collapse of the first star at z = 30. 

These seed black holes then grow in the centers of ~ 
10 6 Mq halo which contains a large amount of high density 
primordial gas, as current theories predict that only one star 
forms per such mini-halo. The dense gas in the central re- 
gion provides abundant fuel for BH accretion. To account 
for their evolution before the major mergers take place, the 
black holes are assumed to grow at the Eddington rate un- 
til their host galaxies enter the simulated merger tree. Such 
an approximation is supported by the fact that the Edding- 
ton ratio in the simulations depends on the galaxy interaction 
and strength of the feedback from the black holes. In our 
simulations, most black holes grow at nearly the Eddington 
rate in the early stages of a galaxy interaction when the feed- 
back is weak. However, when the interaction and the feed- 
back become stronger, the Eddington ratios fluctuate by or- 
ders of magnitude. So a constant accretion rate at the Edding- 
ton limit is no longer appropriate, as we show below. Under 
this assumption, the first progenitor galaxies (Gl and G2) of 



the quasar host have black hole seeds of order 2 x 10 4 /r' M q 
by the time it enters the merger tree at z = 14.4. However, we 
should emphasize that this assumption serves only as an upper 
limit of the early growth of the black holes. Our results in the 
next sections imply that even if all the black hole seeds had a 
uniform mass of ~ 10 5 M Q when they enter the merger tree, 
it is still possible to build a massive one to 10 9 M Q at z ~ 6.5 
through gas-rich mergers. 

In our model, mergers are invoked in the formation of the 
most massive black holes of > 10 7 M Q because that requires 
large supplies of gas. Early on, however, this may not be nec- 
essary to grow the black hole seeds from ~ 100M Q to the 
~ 10 5 M Q we start from, because the accretion rate is small 
so other gas fueling could b e sufficient. As demonstrated in 
Hopkins & Hernquisi] d2006l) . faint AGNs could be fueled by 
stochastic accretion of cold gas that does not involve mergers. 
A similar process could go on in the black hole seeds left by 
the PopIII stars at very high redshifts. We should point out 
that in our simulations, it is necessary for galaxy progenitors 
in the merger tree to have reasonable massive black hole seeds 
(~ 10 5 M Q ) initially in order to build a 10 9 M© black hole at 
z ~ 6.5. However, our results are insensitive to specific forma- 
tion recipes of the seeds. The formation of seed black holes at 
high redshifts is a challenging problem, and some of the pro- 
posed scenarios mentioned above may indeed be necessary to 
make our seeds. However, currently there is no observation 
available to test these models. 

In the picture we adopt in which the seeds come from 
the first stars, the early accretion may be complicated by 
the fee dback from the sta r s. W e note that r ecent stud- 
ies bv IJohnson & Brornml (120061) : lAbel et al.l (120061) and 
Yoshi daet al.l (|2007) show that HII regions form around the 
first stars, and that the halo gas would be photo-ionized, 
photo-heat ed, and evacuated by the radiation feedback from 
the stars. iJohnson & Brom m (2006) suggest that such feed- 
back would deplete the gas in the central region, and would 
delay the black hole accretion by up to 10 8 yrs. However, 
this destruction effect depends sensitively on the lifetime of 
these massive stars, and more importantly on the environment 
which determines both the gas density profile and gas replen- 
ish from inflow of the expell ed gas or neighboring halo s. In 
the simulations presented in IJohnson & Brornml d2006l) . the 
box size is only lOO/r'kpc, too small to contain the large 
scale gravitational potential and the large wavelength density 
modes that drive gas infall, so the initial gas density is low 
and the destruction timescale is long in this case. However, 
the quasar halo in our simulation resides in the highest den- 
sity peak in a volume of 1 h~ l Gpc 3 , where the halo potential 
and gas density, as well we the accretion rate are much higher 
dGao et al.ll2006l) . For a 200 M Q black hole, the accretion rate 
at Eddington limit is only 10~ 6 M Q yr" 1 , which corresponds to 
the Bondi accretion of molecular gas with a typical tempera- 
ture of ~ 100 K at density ~ 10 2 cm" 3 , as implied from equa- 
tion ([8} - Such a density requirement is satisfied with the initial 
conditions of our model. Therefore, the gas re-incorporation 
timescale in our case may be substant ially shorter than that 
estimated in Johnson & Bromm (2006). We will investigate 
in a future project the growth and evolution of the early black 
holes after the death of the first stars in such a cosmological 
environment, using hydro-radiation simulations that include 
both radiative transfer and black hole accretion with ultra- 
high resolutions. 
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2.5.4. Numerical Parameters of Merger Simulations 

The merger tree contains eight galaxies engaging in seven 
major mergers at different times. For each merger event, 
the initial orbits of the incoming progenitors are set to be 
parabolic, consistent with the majority of the major merg- 
ers in our cosmological sim ulation and with previous findings 
dKhochfar & Burkert 2006). The orientation of each merging 
galaxy is selected randomly. The initial separation between 
each merging pair is set to Rq = R w \ r , where R v \ r is the virial ra- 
dius of the incoming system, while the pericentric distance is 
chosen as R p = O.SR^, where R^ is the radial disk scale length 
of the incoming system. We have tested different choices of 
R p and orientations, and found that the impact of these param- 
eters is minor because the orbital properties of the progenitors 
change rapidly through interactions with the multiple galaxies 
in the system. 

Throughout the merger simulation, the mass and force res- 
olutions are fixed for each particle type, and the total initial 
particle number of 1 .0 x 10 results in particle masses of = 
1.1 x 10 7 /i _1 M Q for the halo and m g , s = 2.2 x 10 6 h~ l M Q for 
both the gas and stars. The gravitational softening lengths are 
eh = 60/i _I pc for halo particles and e g s = 30/z _1 pc for both gas 
and stars. In the simulations, it is impossible to resolve indi- 
vidual stars, and the accretion radii of some small black holes 
are under-resolved. However, with the sub-resolution imple- 
mentation in our models, we can calculate time-averaged rates 
of star formation and black hole accretion from the large-scale 
properties of the gas, which are well resolved in our simula- 
tions. Resolution studies of a single merger (Spring e! et al.1 
2005b) with particle numbers from 1.6 x 10 s to 1.28 x 10 7 
show that resolution affects some fine structures of the gas 
and the instantaneous growth rates of star and black holes, but 
the time-averaged properties of the system converge to within 
20%. 

2.5.5. Halo Escape Velocity 

In a galaxy merger with black holes, the black holes may 
merge into one, or may be ejected by gravitational recoil in 
the final stage. Their fate depends on the halo escape ve- 
locity V esc . If the recoil velocity is larger than V esc , then 
the black hole will be k icked out of the halo. We follow 
iBinney & Tremainel(ll987l) to calculate this important param- 
eter y esc . It is defined by 



V esc (r) = V2|*W| : 



(10) 



where <I>(r) is the gravitational potential at a given radius r. 
Because the halo is spherical, the potentials of different spher- 
ical shells add linearly, so $(r) is contributed by two parts, 
i.e., shells within r (/ < r) and outside (r' > r): 



$(r) = -4wG 



i r r°° 

-J p H (rV 2 dr' + j p H (r')r'dr' 



(11) 



where pn(r) is again the iHernquistl d 19901) density profile of 
dark matter halo as in Equation|2] 

Figure [4] shows the escape velocities of the halo progen- 
itors Gl - G8 in Table 1, as well as two merger remnants 
at z — 14 and z — 6.5, respectively. The escape velocity de- 
pends on the halo mass, redshift, and distance from the halo 
center. The V esc remains constant in the central region, be- 
gins to decline around 0AR w \ r . At the center, V esc ~ 2.5V v i r , 
while at the virial radius R vn , the escape velocity is compa- 
rable to the virial velocity (by a factor of ~ 1.5). The iso- 
lated halo progenitors Gl - G8 have a V esc range of ~ 385 
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FIG. 4. — The halo escape velocity V C sc as a function of distance i?/i? v ir 
(/?vir is the virial radius) to the halo center for various models in our merger 
simulations. This plot includes the isolated halo progenitors Gl - G8 in Ta- 
ble 1, as well as the first merger remnant at z — 14 and the last one at z ~ 6.5, 
as labeled in the legend. The shaded region indicates the range of the escape 
velocities of the mergers in our simulations, with the values in the central 
regions being 486kms-' < V esc < 1284kms _1 . 



- 1029 kms" 1 . The first merger halo at z — 14, which has a 
mass of 1.66 x 10 u Mq as the merger of Gl and G2, has a 
central escape velocity of Ve SC ~ 486 kms -1 , while the final 
merger halo at z — 6.5, which has a mass of 7.7 x 1O 12 M , 
has V esc ~ 1284 kms -1 . The shaded region indicates the range 
of the halo escape velocities of the mergers in our simulations. 
In particular, the escape speed in the halo central region has a 
range of 486kms _1 < V esc < 1284kms _1 . This range is im- 
portant for analysis of black hole ejection from gravitational 
recoil in the black hole binaries in § 13. 41 and § 14.21 



3. FORMATION OF A LUMINOUS Z ~ 6 QUASAR 

3.1. Hierarchical Assembly of the Quasar Host 

The vigorous merging history of the quasar host is illus- 
trated through selected snapshots of the gas and stellar dis- 
tributions in Figure [5] and Figure [6j respectively. The pro- 
genitors at high redshifts are very compact and gas rich. 
As the host galaxy of the quasar builds up hierarchically, 
strong gravitational interactions between the merging galaxies 
lead to tidal tails, strong shocks and efficient gas inflow that 
triggers large-scale starbursts, a phenomenon tha t has been 
demo n strated by many numerica l simulations (e.g..|Hernquistl 
[19891: iHernauist & Katzl Jl989l IBarnes & Hernauisll 11991 
1996[lMihos & Hernauis1l994ll996[:ISpringe]||2000l :lBarne 
20021 iNaab & Burked 120031: iLi et al]|2004l) as reviewed by 
IBarnes & Hernq uist ( 1992). The highl y concentrated gas fu- 
els rapid accretion on to the SMBHs dDi Matteo et alj|2005t 
Sprin geTet al.ll2005bh . Between z ~ 14-9, the merging sys- 
tems are physically small and the interactions occur on the 
scale of tens of kiloparsecs. By z ~ 9-7, when the last major 
mergers take place, the scale and strength of interactions have 
increased dramatically. Galaxies are largely disrupted in close 
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FIG. 5. — History of the quasar host shown in selected snapshots. The images give the projected gas density, color-coded by temperature (blue indicates cold 
gas, yellow indicates hot, tenuous gas). The black dots represent black holes. There are eight galaxies in total, engaging in seven major mergers along the timeline 
of the merger events as listed in Table 1. Top panels show interactions in the early stage from z ~ 13 to 9. Middle panels show the last major mergers between 
z ~9— 7, and bottom panels show the final phase. All the galaxies coalesce at z m 6.5, creating an extremely luminous, optically visible quasar (see Figure[6}. 
At this time, there are three black holes, but the luminosity is dominated by the most massive one, which is more than two orders of magnitude larger than the 
others. These black holes merge into a single one at later time. The scale bar indicates a size of 20kpc (comoving), corresponding to an angular size of 3.6" at 
redshift z = 6.5. 
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z- 12.75 z- 10.A2 z- M7 




Fig. 6. — Same as Figure|5] but here the images show the projected stellar density, color-coded by the specific star formation rate (star formation rate per unit 
stellar mass). Blue indicates massive star formation in the galaxies, while red indicates little star formation. To illustrate the quasar activity, we have generated 
"rays" around the quasar. The number and strength of the rays are proportional to the bolometric luminosity of the black holes. These rays are artificial and 
serve only as a visual guide. The systems in top panels are blue, small and perturbed. The quasars appear very faint and buried. In the middle panels, strong 
interactions between galaxies boost star formation and black hole accretion, creating highly irregular morphologies and extremely blue galaxies. The quasars are 
heavily obscured by dense gas. At a later stage (bottom panels), feedback from the black holes quenches star formation, allowing the galaxy color to redden. The 
quasar becomes optically visible as strong outflows blow out the gas. It has a maximum luminosity around z ~ 6.5 when all the galaxies coalesce. After that, 
both the quasar activity and star formation gradually die down, leaving behind an aging stellar spheroid. 
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encounters, tidal tails of gas and stars extend over hundreds of 
kiloparsecs, and intense bursts of star formation are triggered. 

The black holes continue to grow rapidly during this period 
but are heavily obscured by a significant amount of circum- 
nuclear gas. During galaxy mergers, the black holes follow 
their hosts to the center of the system and can interact closely 
with each other. It has been shown that black hole binaries 
decay rapi dly in a gaseous env iro nment and can merge within 
~ 10 7 yrs dEscala et al.ll2004t O I20071) . Because the galax- 
ies in our simulations are very gas rich and the gas is highly 
concentrated during the mergers, we therefore assume that the 
black holes merge efficiently owing to st rong dynamical fric- 
tion with the gas (ISpringel et alJl2005bl) . We wil l return to 
more discussions of this process in § !3.4l and § 14.21 

At redshift z ~ 6.5 the progenitor galaxies coalesce, induc- 
ing high central gas densities that bring the SMBH accretion 
and feedback to a climax. The SMBH feedback then drives 
a powerful galactic wind that clears the obscuring material 
from the center of the system. The largest SMBH b ecomes 
visible as an optically-bright quasar (Hop kins et al.1 [2006a) 
during this phase, after which quasar feedback quenches star 
formation and self-regulates SMBH accretion. Consequently, 
both star formation and quasar activity die down, leaving a 
remnant which reddens rapidly, as illustrated schematically in 
Figure [6] 

3.2. Star Formation History 

The evolution of the star formation rates (SFRs) of each 
individual galaxy, and total SFR of the whole system are 
shown in Figure|7](fo/9 panel). The system forms stars rapidly 
as these compact and gas rich progenitors undergo strong 
interactions. 
> 10 4 Moyr 

begin their final major mergers, while the SFRs of individual 
galaxies fall below a fewx 10 3 M Q yr" 1 , within the starburst 
intensi ty li mit of 10 3 M^, yr" 1 kpc~ 2 proposed by Meur er et al.l 
d 19971) and lThompson et al.1 (|2005). At z < 7 the star forma- 
tion rate decreases gradually owing to a depletion of the gas 
supply and progressively stronger feedback from the SMBHs. 
At the time of final coalescence (z ~ 6.5) the star formation 
rate is ~ lOOMQyr -1 , an order of magnitude lower than for 
estim ates of Jl 148+5251 dBertoldi et al.ll2003at ICarilli et al.l 
2004). We note, however, that the estimates by these authors 
are based on the assumption that the FIR luminosity is domi- 
nated by young stars, and they cannot rule out the possibility 
that AGN may contribute sig nificantly to th e luminosity. 

In a forthcoming paper (ILietal.l[2007h . we have calcu- 
lated the infrared properties of the quasar system using a 3-D 
Monte Carlo radiative transfer code that incorporates adap- 
tive grids and treats dust emission self-consistently. We find 
that the far-infrared luminosity of our quasar is not domi- 
nated by young stars but instead has a substantial quasar con- 
tribution of over 80%. This finding is supported by obser- 
vations of Jl 148+5251 i n near-IR (e.g., Charmandar is et all 
12004 iHines et al.l 120061 and more recently iDwekl 120061) . 
which show a remarkably flat spectral energy distribution and 
suggest an AGN origin for the flux exc ess. Furthermore , 
adopting a total gas ma ss of ~ 10 I0 M© (IWalter et al J 120041: 
iNaravanan et al.ll2006d) in Jl 148+5251, a si mple application 
of the Schmidt- Kennicutt star formation law (Schmidt 1959; 
lKennicuttl [l998 ) gives a star formation rate of ~ 2OOM0vr -1 , 
close to what we find here. 

Within only about 600 Myrs from z = 14.4 to z = 6.5, the 
system accumulates a stellar mass of ~ 10 12 M Q as shown 



The total SFR ranges from — lOOM yr _1 to 
between redshifts z ~ 9-8 when the galaxies 
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FIG. 7. — Time evolution of star formation rate (SFR, top panel) and spe- 
cific star formation rate (SSFR, bottom panel), respectively. The colored lines 
indicate individual galaxies, while the black lines give summed quantities for 
the entire system. 



in Figure [7] (bottom panel). The specific star formation rate 
(SSFR), or "b-parameter" (e.g., iBrinchmann et al.ll2004l) . is 
defined as SSFR = SFR/M star . It is a measure of the fraction 
of the total stellar mass currently forming at a specific time. 
The SSFR is typically larger in high-redshift galaxies than in 
ones at low redshifts owing to vigorous star formation. Dur- 
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FIG. 8. — Time evolution of mass-weighted metallicity in the quasar host, 
from gas (blue curve), stars (red curve), and the mean value in the central 
region (R < lkpc) of each galaxy (black curve). The blue hatched region 
indicates the range of 25%-75% of the gas metallicity. 



ing the past several years, there has been rapid progress in ob- 
serving galaxies at z > 6 using the Hubble Space Telescope 
(HST), and the Spitzer Space Telesc ope (Spitzer) coupled 
with ground-based observatories (e.g ., [Dickinson et al ]l2004t 
Bunke retalJ I2004I: iBouwens et alJ l2004t Giavalisc o et al. | 
20041: lEgami et all I2005I: lEyles et all I2005I: iMobasher et all 
200llYan et al.H2005L I2006t Evles et al.ll2006l) . and hundreds 
of these distant objects have been detected. These frontier 
observations suggest that the Universe experienced rapid star 
formation during the redshift interval 14 > z > 6, and the de- 
velopment of large stellar systems in th e mass range of ~ 
10 10 ~ 11 Mp. In particular, several groups (E2armetaLlf2005; 
lYan et al.ll2006tlEvles et a"T1l2006l) find SSFRs in the range of 
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FIG. 9. — Spatial distribution of mass-weighted metallicity of the quasar 
host at z ~ 6.5, from both gas (top panel) and stars (bottom panel), respec- 
tively. The images are projected metallicity adaptively smoothed over 32 
particles (analog to the SPH kernel in a 2-dimensional plane). The black dot 
indicates the center of the quasar. 



3.3. Metal Enrichment 

Rapid star formation in the quasar progenitors produces an 
abundant mass of heavy elements to enrich the ISM. Obser- 
yations of Jl 148+525 1 show solar metallicity in the system 
dBarth et alJ [200l IWalter et alJ 12001 iMaiohno et all l200l 
iBecker et al.l2006l) . Figure[8]shows the metallicity in our sim- 
ulated quasar system at different times. Note that the dips 
and jumps in the curves owe to incoming new galaxies which 
bring in metal-poor pristine gas and newly-formed stars. The 
quasar host reaches solar metallicity as early as z ~ 12, and 
maintains similar levels to later times. The spatial distribution 
of metallicity from both gas and stars at the peak quasar phase 
at z ~ 6.5 is shown in Figure[9] The metals are widely spread 
owing to outflow from the quasar feedback and gas infall to- 
ward the merger center. The metallicity in the central region 
of the merger remnant is slightly above the solar value. In 
some outer regions, because the gas and stars are still falling 
back to the system center, the infalling material triggers small- 
scale bursts of star formation. So the metallicity in these blobs 



appears to be super-solar, as shown in Figure[9] 

Calculations of carbon monoxide emission using non- 
local thermodynamic equ ili brium radiative t r ansfer codes 
dNaravanan et al . 2006a b) by Narayan an et alJ (l2006cl) show 
CO luminosities, excitation patterns, and morphologies 
within the central ~ 2kpc of the quasar host center that 
are consistent with observa t ions of Jl 148+525 1 (IWalter et al.1 
l200llBertoldTetd]l2003bl: IWalter et alJl2004h . These results 
were derived using Galactic CO abundances, and thus support 
our conclusions that significant metal enrichment takes place 
early in the quasar host, as a result of strong star formation in 
the progenitors. 

3.4. Growth of Supermassive Black Holes 

In the simulations, the quasar host at z ~ 6 is built up by 
eight progenitors, each containing a black hole in the center. 
Figure[l0] shows the evolution of the black hole accretion rate, 
the Eddington ratio, and the integrated masses of the whole 
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FIG. 10. — The growth history of the quasar system, including the black 
hole accretion rate BHAR, top panel, the Eddington ratio Lbol/^Edd (middle 
Panel), and black hole mass (bottom panel). Note that the black curves rep- 
resent totals, while colored curves show individual black holes, as indicated 
in the legend. 



system and individual black holes. The total black hole ac- 
cretion rate grows steadily during the hierarchical assembly 
of the host galaxy and peaks at ~ lOM yr -1 around z ~ 6.5 
during the final coalescence. 

The Eddington ratio, Lboi/^Edd, of each individual black 
hole varies with time, depending on the galaxy interaction and 
feedback from the black holes. The black holes maintain ac- 
cretion at the Eddington limit for only a fraction (< 50%) of 
the time. At the peak of quasar activity, the Eddington ra- 
tio of the most massive black hole is near unity, while that of 
the other black hole is only 0.1. However, collectively, the 
whole system appears to accrete at Lboi/^Edd ~ 1 at z > 6.5, 
as implied in Figu r e ITOl Studies of black hole accretion (e.g., 
Vestergaardll2004t iKollmeier etai1l2005t IVestergaard & Fan! 
2006) show that the Eddington ratio has a wide range of 
0.01-1.0, and it varies with both luminosity and redshift. 
Luminous systems tend to have higher Lboi/^Edd than less- 
luminous counterparts, and at z > 4, most quasars shine at 
nearly Eddington luminosity. Our results suggest that individ- 
ual black holes do not always necessarily accrete at the Ed- 



dington rate. However, since high-redshift, luminous quasars 
may form through mergers of several galaxy progenitors con- 
taining black holes as in our case, therefore the growth of 
the quasar represents a collective contribution from each in- 
dividual black hole. The total black hole mass increases from 
- 6 x 10 4 M Q at z « 14 to about 2 x 10 9 M o at z » 6.5, c lose 
to that estimated f or Jl 148+5251 bv lWillott et ail d2003l) and 
iBarthet al.1 (120031) . 

In the simulations, we do not have sufficient resolution nor 
the relativistic physics to consider the ejection of black holes 
by gravitational recoil. The black holes are assumed to merge 
efficiently once their separation is below the spatial resolu- 
tion. In the final stage of black hole mergers, the emission 
of gravitational wave carries linear mo mentum, which could 
cause the black h oles to recoil (e.g., iBonnor & Rotenbergl 
I196U lPeresl[l962r) . If the recoil velocity is larger than the 
halo escape velocity, the n the black holes will be kicked 
out from their halo (e.g.. iFitchettl [l983h iFayata et all 1 2004 ; 
Merritt et al.ll2004t IMadau & Ouataert 2004: IHaimanl 12004: 



Yoo & Miralda-Escude 2004; Volonteri & Rees 2005). Previ- 



ous studies by (|Haimanll2004t lYoo & Miralda-Escude] 120041 : 
IVolonteri & Reesll2005t IHaimanl 120061) suggest that constant 
or super-Eddington accretion is required to produce 10 9 M Q 
black holes a t z ~ 6 if eject ion of black holes is included. 
In particular, Haiman (2004) suggests a black hole will be 
ejected if the kick velocity V^ick for the coalescing SMBH bi- 
nary is larger than twice the halo velocity dispersion Chaio, 



ick 



> 



2CThaio, as the dynamical friction timescale for the 
kicked black hole to return to the halo cent er is longer than 
the Hubble time ( IMadau & Ouataertl 12004). By applying 
this ejection criterion to a Press-Schechter merger tree of a 
8.5 x 1 12 Mfp halo w ithin which Jl 148+5251 is assumed to 
reside, IHaimanl d2004l) finds that the SMBH of the quasar 
gains most of its mass rapidly from seed holes during 17 < 
z < 18 due to black hole ejection, and the SMBH likely ac- 
cretes with super-Eddington rate in order to build a mass as 
that of J 1148+5251. 

However, the halo escape velocities or velocity dispersions 
in our model are much larger than the currently best esti- 
mates of the kick velocity. The quasar halo in our simu- 
lations has active merging history from redshifts z — 14.4 
to z — 6.5, and the halo progenitors have mass es much 
higher than those considered in the previous studies dHaimanl 
120041: lYoo & Miralda-Escud^l2004t IVolonteri & ReeslHooH) . 
The quasar halo builds its mass from ~ 1.16 x 1O U M0 at 
z ~ 14 (the sum of progenitors Gl and G2 in Table 1) to 
7.7 x 10 12 M Q at z = 6.5. As shown in Figure [4] (the shaded 
region), the central escape velocity of the mergers in our 
simulations is in the range 486 - 1284 kms -1 . Currently, 
the maximum kick velocity for unequal-mass, non-rotating 
BH binary is in the range of ~ 74 - 250k ms~ from both 
analytic post-Newtonian approxima tion (e.g., iBlanchet et alj 
120051: [bamour & Gopakumar 2006) and the ground-breakin 
full r e lativistic numerical simulations (e ,g.. | Herrmann et al 
120061: iBakeretalJ [20061: IGonzalez et al]l2006l) For eaual- 
mass, spinning BH binary. iFavata et alj (120041) estimate a 
ran ge of ~ 100 - 200km s " 1 using BH perturbation theory, 
and iHerrmann e t al. (2007) derive a formula from relativis- 
tic simulations, Vyck = 4755km s~ , where S < 1 is the BH 
spin. This gives a maximum kick of 475kms~ for max- 
imal spin. Although it is also reported that the recoil ve- 
locity can be as large as tho usands kms -1 (IGonzalez et alj 



ty can be as large as th ousands kms ouonzalez et alj 
3 ICampanelli et all 120071) for BH binary in the orbital 
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FIG. 11. — Evolution of total BH (black curve) and stellar (yellow curve) 
mass, respectively. The stellar mass is multiplied by a factor of 0.002, to 
reflect the observe d correlation of Mbh ~ 0.002M stal at the present day, as 
parameterized by Marconi &Hunl|j2003l) . 



pla ne with opposite-directed spin. However, as pointed out 
by iBogdanovic et alj d2007l) . such a configuration is rather 
uncommon, especially in gas-rich galaxy mergers, because 
torques from accreting gas suffice to align the orbit and spins 
of both black holes with the large-scale gas flow. The re- 
sulting maximum kick velocity from such a configuration is 
< 200 km s -1 . Overall, the kick velocity from the latest cal- 
culations of black hole binary is in the range of ~ 100 - 
475km s -1 , falling safely below the escape velocities of the 
quasar halos in our simulations, so black hole ejection may be 
insignificant in our case. 

Moreover, we find that our model can produce a 10 9 M Q 
SMBH even if ejection is allowed. From Figure [TOl (bottom 
panel), the 10 9 M© SMBH is dominated by BH5, most of its 
mass comes from gas accretion. Even if the less massive black 
holes, for example BH7 or BH8 was ejected, the most mas- 
sive one BH5 is still able to reach ~ 10 9 M Q in the end. Fur- 
thermore, even if all the seeds started with ~ 10 5 M Q in the 
merger tree, the result would be about the same. We therefore 
conclude that the results from our modeling are robust. Sup- 
permassive black holes of ~ 10 9 M Q can grow rapidly through 
gas accretion and mergers hierarchically in the early Universe, 
constant or super-Eddington accretion is not necessary, unless 
the recoil velocity of the coalescing black hole binary is ex- 
tremely high such that most of the black hole seeds in our 
simulations are ejected (e.g., Vki C k > lOOOkms -1 ). 

3.5. Correlations between Supermassive Black Hole and 
Host Galaxy 

Tight correlations between supermassive black 
holes and hosts have been o bserved in local galaxies 
(e.g., iMagorrian et alJ 119981: iFerrarese & Merritf 1200(1 
iGebhardtet al. 2000), but the inference of these relationships 
at higher redshifts remains an open question. Because the 
eight galaxies in our system interact vigorously with each 



other, the stellar components are widely spread and mixed, 
it is impossible to separate individual galaxy-SMBH pairs, 
so we only consider the correlations in total quantity of 
the whole system. A comparison of the total stellar mass 
and total black hole mass is shown in Figure QT| At early 
time, both the stars and black holes grow rapidly through 
galaxy mergers. Shortly after the peak quasar phase, strong 
feedback suppresses both the accretion and star formation, 
the masses of the black holes and stars become saturated 
gradually, and in the end satisfy Mbh ~ 0.00 2M star , similar to 
the correlation measured in ne arby galaxies (Ma gorrian et alj 
119981: iMarconi & Hunt! 120031). Our res ults are consistent 
with fi ndings by Robertson et al. (2006b]) and iHopkins et alJ 
d2007l) . and demonstrate that the observed M B H-Mbui g e 
scaling relation is a result of the coeval growth of the SMBH 
and its host galaxy, and that it holds across different cosmic 
times. 

We note, however, that the velocity dispersion of the 
stars in the remnant center is about ~ 500kms _1 (after 
the system relaxes) owing to the deep potential well of 
the merger system, so the Mr h-ct relation falls below 
the correlation observed locally (IFerrarese & Merrittl 120001 : 
iGebhardt et alj 120001 : iTremaine et alj 120021) . Single merg- 
ers of pro genitor galaxies construc ted in a redshift range of 
z=0-6 by Robertson et alj d200 6b) appear to follow the ob- 
served Mbh-c correlation with a weak redshift dependence 
of the normalization, which results from an increasing ve- 
locity dispersion of the progenitors at higher redshift. The 
multiple mergers we derive from cosmological simulations 
take place at much higher redshifts and hence the progeni- 
tors have larger velocity dispersions, implying a larger de- 
viation from the local M bh~Q" relation than in the work of 
iRobertson et al.l d2006bl) . However, because we do not follow 
subsequent mergers and accretion into the host halo below 
z ~ 6, the implications of this result for the evolution of the 
Mbh-c relation are unclear. 

Observations of active galaxies have yielded ambiguous re- 
sults about th e SMBH-spheroid relationship. For example, 
iGreene & Hoi ((2006) report a lower zero-point of the Mbh-c 
of local active galax ies than t hat of the inactive sample 
dTremaine et aUl2002h : at z > 0, IShields et al.1 d2003h found 
the same M rh~Q" relation in t he redshift range z =1-3, while 



others (e.g.. iTreu et al.l 
l200llPengetal]l2006t 



20041 IWalter et al 



12004 iBorvs et alj 



Shields et al. 2006) show correlations 



with various offsets. In particular, IWalter et alJ (l2004h es- 
timate a dynamical mass of ~ 5 x 10 10 Mq using the CO 
linewidth measured in J 1 148+525 1 , and suggest that the bulge 
is under-massive by at least one order of magnitude compared 
to the l ocal MBH~Mhui gP relation . However, the CO calcula- 
tion by Naravanan et al. (2006c) finds that the CO linewidth 
of the quasar in our simulation is l arger th an t he mean 280 
km s" 1 measured by iBertoldi et al.l d2003bl) and IWalter et alJ 
(2004) by almost an order of magnitude, and that the derived 
dynamical mass is ~ 1Q 12 M^, putting the simulated q uasar 
on the MBH _ Mb u ige correlation. Narayanan et al. (2006c) fur- 
ther suggest that the observed emission line may be sitting on 
top of a much broader line, which may be tested by future 
observations with large bandwidths. 

The different relations reported from the observations may 
reflect a divergence of the methods used to estimate the black 
hole mass and stellar properties, o r may represent different 
evolu tionary stages of the systems ( Wu 2006; Hopki ns"et alJ 
120071) . More observations and better measurements of black 
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hole mass and properties of host bulges will be crucial to 
study the SMB H-host relations in high-redshift quasar sys- 
tems ( Vestergaard & Peterson 2006) and to test our hypothe- 



3.6. Quasar Luminosities 

Both the bolometric and attenuated luminosities of the 
quasar and the host galaxy in the simulation s can be read- 
ily calc ulated following the methodology of Hopki nTet al.l 
(l2005dl) . The bolometric luminosity Lboi of stars is cal- 
culated using the stella r population synthesis model of 
iBruzual & Charlotl d2003l) . while that of a black hole is cal- 
culated as Lboi = £ r Mc 2 , where e r = 0.1 is the radiative effi- 
ciency, M is the black hole accretion rate, and c is the speed 
of light. In this calculation, the black holes are assumed to 
be non-rotating. If the black holes are spinning, their radia- 
tive efficiencies and luminosities would be higher due to the 
shrink of the innermost stable circular orbit, by up to a factor 
of 4 for maximal rotation. 

The B-band luminosity of each source is corrected for at- 
tenuation by absorption from the ISM along the line-of-sight. 
We first calculate the line-of-sight column-density of the gas 
from each source to a distant observer. For each black hole 
we generate 1000 radial sight lines originating at the black 
hole particle location and uniformly spaced in the An solid 
angle dcos0d<f>, while for the stars, an accurate estimate of the 
luminosity is possible with only one sight line per source ow- 
ing to the extended distribution. Along each ray, the gas col- 
umn density is calculated using a radial spacing of Ar = T]h sm \, 
where 77 < 1 and h sm \ is the local SPH smoothing length. The 
distribution of line-of-sight properties converges for > 100 
rays and at a distance of > lOOkpc. In the calculation, only 
the diffuse-phase density is considered because of its large 
volume filling factor > 99%, allowing for a determination of 
the lower limit on the column density along a particular line 
of sight. 

Adopting the mean obse rved intrinsic quasar continuum 
spectral energy distribution (Richa rds et al J 12006) gives a B- 
band luminosity which is well approxim ated by the following 
equation given bv lMarconi et all (120041) . log (L bo i/L B ) = 0.80- 



0.067L + 0.017L 2 -0.0023L 3 , where L = log(L bo i/L Q ) - 12, 
and Xb = 4400 . We then use the Milky Way gas-to-dust ratio 
scaled by metallicity, A B /N H = (Z/0.02)(A b /Nh)mw to deter- 
mine the extinction along a given line of sight for this band. 
In the above calculation, we do not include a full treatment of 
radiative transfer, and therefore do not model scattering or re- 
processing of radiation by dust in the infrared. However, for 
the B-band luminosity, results using a 3-D Monte Carlo radia- 
tive transfer code are close to those calculated using the met h- 
ods we present here dLi et al.ll2007t IChakrabarti et alj |2006). 

Figure[l2]shows both the bolometric and attenuated B-band 
luminosities of the quasar, compared with observations of 
J 1 148+525 1 . The system is intrinsically bright with a total lu- 
minosity > 1O U L0, and the host appears as an ultraluminous 
infrared galaxy (ULIRG) with Lboi > 10 12 L© for most of the 
time. At high redshifts, z > 8, starlight dominates the total lu- 
minosity. However, black holes take over at a later time. The 
quasar light-curve increases dramatically, peaking at z ~ 6.5, 
when it is powered by the most massive black hole accreting 
at near the Eddington rate. The estimated Lboi of J 1 148+5251 
differs from that of the simulated quasar by less than the un- 
certainty in the luminosity estimate. The rest-frame B-band 
absolute magnitude reaches Mb ~ -26.5, almost one mag- 
nitude fainter than that of J 1148+5251 derived from 1450 



data dFan et al.ll2003l) . However, we should emphasize that 
in this paper, our main goal is to investigate the plausibility of 
forming luminous z ~ 6 quasars through hierarchical mergers, 
rather than precisely reproducing the properties of an individ- 
ual quasar such as Jl 148+5251, so the disagreement shown 
in Figure [T2] should not be taken too literally. Moreover, the 
exact luminosity can change by a factor of several from rela- 
tively trivial or random details in the simulations. If the black 
hole spin is taken into account, then the simulated luminosi- 
ties would increase by a factor of up to 4, which would match 
the observation of J 1 148+5251 better. 

Feedback-driven outflows create un-obscured lines-of- 
sight, allowing the growing central SMBH to be visible as an 
optically-bright quasar between redshifts z ~ 7.5 and z ~ 6.4. 
At the peak of the quasar activity, more than 50% of the 
1000 sight lines have Lb > 10 12 L Q . The absorbed light is 
re-emitted at infrared wavelength s by dust. W e find that 
the luminosity in the far infrared (Li et al. 2007) is close to 
L F i R ~ 10 13 L W estimated for Jl 148+5251 by iBertoldi et all 
d2003al) . Moreover, we find that up to 80% of the FIR light 
comes from the black hole, while stars contribute only ~ 20%. 
This may explain why the star formation rate at z w 6.5 during 
the peak quasar phase is an orde r of magnitude lower than the 
estimate from FIR observations (Bertoldi et al. 2003a), which 
will be contaminated by the AGN. 

Another prominent feature of Figure [12] is a clear 
phase transition from starburst to quasar. It has long 
been suggested that ULIRGs are pow ered by starbursts 
in galaxy merger s (for reviews, see ISanders & Mirabel 
1996; Joged 120061) , and that bright quasars are the descen- 
dants of ULIRGs (ISanders et al.l 119881; iNorman & Scovilld 
19881 IScovilld 120031: lAlexander et al.l 120051) . This" conjec- 
ture has been supported by observations of quasar hosts 
(e.g. IStocktonl 19781; iHeckman et al.lll984tlHutchings & Neffl 
1992; Bah calletal.lll997UHutchingsll2005l). and theoretical 
modeling (Ho pkins et alJl2006al) . In iLi et afl (120071) . we cal- 
culate the spectral energy distributions (SEDs) of the quasar 
system and its galaxy progenitors. We find that the SEDs of 
the syste m at z > 8 are character ized by those of starburst 
galaxies (Sanders & Mirabel 1996), while at the peak quasar 
phase, the SEDs resemble those observed in z ~ 6 quasars 
(Jiang et al. 2006). We also find that the system evolves from 
cold to warm ULIRG as it transforms from starburst to quasar 
phase. Our results provide further theoretical evidence for 
the ULIRG-quasar connection in quasar systems in the early 
Universe. 

The quasar lifetimes d epend on the observe d luminosity 
threshold, as proposed by iHopkins et al.l (l200 5d). In our sim- 
ulation, at the peak luminosity of Lb w 2 x 10 12 L Q , the quasar 
lifetime is roughly ~2x 10 6 yrs, as shown in Figure [13] 
Again, If black hole spin is included in the calculation, the 
luminosity of the quasar would increase by a factor of sev- 
eral, and the quasar lifetime would be longer. However, when 
increasing the radiative efficiency, the Salpeter time (e-folding 
time for Eddington-limited black hole growth, Salpeter 1964) 
is increased by an identical factor, meaning it would also re- 
quire a longer time to reach the same mass. If high-redshift 
quasars are rapidly rotating, then, our calculations demand 
that either the seeds be much more massive at z > 6, or that 
they accrete in a super-Eddington manner. In other words, if 
the observed Sloan quasars at z ~ 6 shine with Eddington lu- 
minosity but are rotating rapidly, then our model suggests that 
their masses would be considerably smaller than e stimated. 

We note that recent Spitzer observations by Jiang et al. 
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FIG. 12. — Comparison of luminosities from simulations and observations. Shown are bolometric (top panel) and attenuated luminosities in the rest-frame 
B-band (bottom panel). Note Lbol of SDSS Jl 148+5251 (green dot) is an estimate for a SMBH of 3 X 10 9 Mq accreting at the Eddington rate, with the error 
bar indicating the mass range of 1-5 X 10 9 Mq I Willott et al. 2003; Barth et al. 2003), while the Lb is converted from observations at wavelength 1450 jFan et al] 
2003). The yellow, red, and black curves represent luminosities of stars, black holes, and total (sum of the above two), respectively. For the black holes, LBH.mean 
is the average luminosity over 1000 sight lines. Note in the luminosity calculation, the black holes are assumed to be non-rotating. If the black holes are rotating, 
their radiative luminosities could be higher by up to a factor of a few, see text for more discussions. 
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FIG. 13. — Quasar lifetimes as functions of different B-band limiting lu- 
minosities. The black and red curves represent the total luminosity of the 
system, and mean luminosity over 1000 sight lines of the black hole, respec- 
tively, corr esponding to the curves of the same color in the bottom panel of 
Figure [T2] 



(2006) show that 2 out of 13 quasars at z ~ 6 have a remark- 
ably low NIR-to-optical flux ratio compared to other quasars 
at different redshifts, and these authors suggest that the two 
quasars may have different dust properties. According to our 
model, however, these two outliers may be young quasars that 
have just experienced their first major starburst but have not 
yet reached peak quasar activity, so the light from star forma- 
tion may be dominant, or comparable to that from the accret- 
ing SMBH still buried in dense gas. This may explain the low 
NIR flux, as well as the B-band luminosity and the narrow L a 
emission line, which are primarily produced by the starburst. 
We will address this question further in a future pap er with 
detailed modeling and IR calculations dLi et a l. 2007). 

4. DISCUSSION 

4.1. Comparison with Previous Models and Robustness of 
Our Results 

Our multi-scale simulations that include large-scale cos- 
mological N-body calculations and hydrodynamic simula- 
tions of galaxy mergers, as well as a self-regulated model 
for black hole growth, have successfully produced a lumi- 
nous quasar at z ~ 6.5 with a black hole mass of ~ 2 x 
10 9 /t'Mq and a number of properties similar to those of 
Jl 148+5251, th e most distant quasar detected at z = 6.42 
dFan et al.l 12003). O ur approach differs from previous semi - 
analytic studies by Haiman & Loebl (12001b: iHaimanl d2004l) : 
Yoo & Miralda- Escude (2004): IVolonteri & Reesl d2005l) and 
Volonteri & Reesl ([2006 ) in the following ways: 



1. We use a realistic merger tree derived directly from 
multi-scale, high-resolution cosmological simulations. 
The previous studies used merger trees of dark mat- 
ter halo s generated with the extended Press-Schechte r 
theory dPress & Schechteill 19741 lLacev & Colelll993l) . 



which may underestimate the abundance of high-mass 
h alos by up to one order of magnitude, as shown in 
§ 12.31 Also, the merger trees in those studies started 
from much higher redshifts than what we consider here. 
In our model, the quasar halo is the largest one in a vol- 
ume of l/T 3 Gpc 3 . It has a mass of ~ 7.7 x 10 12 M© 
at z ~ 6.5 built up through seven major mergers from 
z~ 14.4 to z~ 6.5. 

We follow the evolution of the system and treat 
the gas dynamics, star formation, and black hole 
growth properly. This approach is critical to inves- 
tigation of the properties of both b lack holes and 
host galaxies, and the i r evolu ti on (e.g.,[ Di Matte o et al.l 
2005; Springel 'etai] l2005ri iRobertson et all 12006b; 
Hopkins et al. 2006a), but it was not included in those 
previous studies on formation of z ~ 6 quasars. 

3. We employ a self-regulated model for the growth of su- 
permassive black holes, in which the accretion is regu- 
lated by the black hole feedback, and the rate is under 
the Eddington limit. In the previous studies, the black 
hole growth was unregulated, but instead a constant- or 
super-Eddington accretion rate was used. 

4. In our simulations, we do not consider black 
hole ejection caused by gravitational recoil ow- 
ing to insufficient resolution and lack of relativis- 
tic physics. However, the halo escape veloci- 
ties in our simulations are in the range of 486 - 
1284kms~ , much larger than the kick velocity ~ 100 
- 475 km s" 1 (e.g.. iHerrmann et al.l 120061: iBaker et alj 
2 006; G o nzalez etal.ll2006tlHerrmannet al.ll2007l see 
§ 12.5.51 and § 13.41 for more details). Therefore, black 
hole ejection may be negligible in our case. Previous 
studies had much smaller halo progenitors at higher 
redshifts than ours, so the black hole seeds would be 
more likely subject to ejection from their halos. This 
leads to the conclusion in these studies that constant or 
super-Eddington accretion is needed owing to signifi- 
cant black hole ejection. Our results are robust within 
the best estimates currently available for the recoil ve- 
locity of the black hole binary. 

5. The black hole seeds in the galaxy progenitors in our 
simulations are massive (e.g., ~ 10 5 M Q at z ~ 14). The 
sub-resolution recipe in our model does not allow us 
to resolve the actual formation and accretion of such 
black holes below this mass scale. The formation of 
these seeds is an unsolved problem, but our results do 
not depend on the specific prescription of the formation 
process. We adopt a picture in which the seed holes 
come from the remnants of the first stars (which have 
a mass 200M Q at z = 30) and grow under Eddington 
limit until they enter the merger tree we simulated. If 
the growth is de layed by radiation feed back from the 
PopIII stars (e.g.. lJohnson & Bromml2006l) . then super- 
Eddington accretion, or other proposed scena rios (e.g., 
lBromm& Loebl 120031: iBegelman etafl 120061) may be 
necessary to form massive seeds of ~ 10 5 M Q in the 
protogalaxies. 

Overall, we conclude that the results from our simulations, 
which are more realistic and more detailed than the mod- 
els previously done, are robust. Suppermassive black holes 
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of - 1O 9 M can form rapidly through gas-rich hierarchical 
mergers under Eddington limit, even within a short period of 
time. We find that constant or super-Eddington accretion is 
not necessary unless the above assumptions in our modeling 
break, i.e., there are no massive black hole seeds of 10 5 M Q 
available at z ~ 14, or the recoil velocity of the coalescing 
black hole binary is extremely high (e.g., > lOOOkms -1 ). 
Under these extreme circumstances, some "exotic" processes 
such as super-Eddington accretion may be necessary to grow 
a ~ 10 9 Mq SMBH within a few hundred mil lion years. How- 
ever, w e should note, as pointed out by iBogdanovic et al.l 
(2007), that most gas-rich galaxy mergers have a configu- 
ration such that the orbit and spins of both black holes are 
aligned with the large-scale gas flow owing to torques from 
accreting gas. Such a configuration has a maximum kick ve- 
locity < 200 kms -1 , which is well below the escape velocity 
of a 10 10 M Q dwarf galaxy, as well as those of the halos in our 
modeling. 

4.2. Merging History of Black Holes 

During the galaxy mergers, the black holes follow their host 
halos to the system center and can form binaries (or mul- 
tiple systems). The coalescence of a black hole binary in- 
clud es three distinct phases: in spiral, merger, and ringdown 
(e.g jFlanagan & Hu ghes 1998). Whether black hole binaries 
can coalesce on short timescales is a matter of debate. In a 
stellar environment, it has been argued that a binary hardens 
very slowly owing to an eventual depletion o f stars that cause 
the b i nary to lose angular momentu m (e.g., Begelma n et al.l 
1980; Milosavljevic & Merritt 2003). In a gaseous environ- 
ment, h oweve r, numerical simulations bv lEscala et al.l d2004l) 
and O d2007l) show that the binaries decay rapidly owing to 
strong dynamical friction with the gas, and they likely merge 
within 10 7 years. Because our galaxies are very gas rich 
and have large central concentrations of gas during the merg- 
ers, we assume the black hole particles coalesce once their 
separation decreases below our spatial resolution (30/i~'pc) 
and their relative spee d falls below the local gas sound speed 
dSpringel et al.ll2005bl) . 

In the simulations, we do not have sufficient resolution nor 
the relativistic physics to consider the ejection of black holes 
by gravitational recoil during the merger phase. However, as 
discussed in § 12.5.51 and § 13.41 the halo escape velocities in 
our simulations are much larger than the maximum kick ve- 
locity for black hole binary estimated from the latest relativis- 
tic calculations. So black hole ejection is likely unimportant 
in our modeling. To accurately address gravitational recoil 
in the galaxy merger simulations, we need to include general 
relativity, resolve the dynamics of black hole binaries with 
extremely high resolution, and calculate the halo potential 
in a cosmological context (in which halo potential distribu- 
tion may be different from that of a single object). However, 
such a comprehensive treatment is impossible at the moment. 
We therefore assume the black holes merge quickly once they 
reach the stage of gravitational radiation. 

These coalescing supermassive black holes will be strong 
sources of gravitational radiation de tectable by the Laser In- 
terferometer Space Antenna (L ISA. lFolknerlll998l). as sug- 
gested by many authors (e.g., [Thorne & Braginskii 1 19761 : 



Haehnelt y i998; Flanagan & Hashes 1998; Menou el al. 2001: 



Hughes! 120021: ISesana et alJ 12005b iKoushiappas & Zentnerl 
2006). By tracing the merging history of the SMBHs, LISA 
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FIG. 14. — Evolution of star formation and black hole growth in merger 
simulations with a starburst-driven wind model. The simulation is run with 
lower resolution (Ntot ~5x 10 5 ), and the specifications of the wind model 
is: wind efficiency r\ = 0.5, wind energy coefficient \ = 0.25, wind free travel 
length L w = 20kpc, and a wind velocity V w = 418kms~' . 



of the associated dark matter halos. Because luminous, high- 
redshift quasars are likely sites of vigorous hierarchical merg- 
ers, they may be the best targets for LISA to explore the early 
Universe. 

4.3. Feedback from Starburst-driven Winds 

Vigorous star formation would induce a galactic wind and 
mass outflow, a phenomenon that has been observed to prevail 
in both local star-forming gala xies as indicated by blueshifted 
optica l absorption l i nes (e.g.. [Martini 1 1999]: | Heckman et alJ 
I2000t iMartinl 120051 iRupke et all I2002L 120051) . and Lyman 
Break Galaxies at z ~ 3 as indicated by blueshifted interstel- 
lar absorption lines and redshifted Ly q emission lines (e.g., 



iPettini et aLll20q2t IShaplev ■ et alJ 
ters at z ~ 5.7 (lAjiki et al.ll2002l) . 



2003), as well as Lya emit- 
These galactic winds are 



could shed light on the distribution, structures and evolution 



generally thought to play a sig nifi cant role in galaxy evolution 
(e.g., see Veille uxet al.ll2005l and lCox et al.ll2006al for recent 
reviews). 

The strong starburst preceding the major quasar phase in 
our simulations may drive strong galactic winds and affect the 
black hole growth. To investigate the impact of the feedback 
from a starburst-driven wind on the growth of the black hole, 
we have done the same merger simulation with lower resolu- 
tion (Ntot ~ 5 x 10 5 ) and with a canonical wind model from 
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FIG. 15. — Comparison of halo abundances at z ~ 6 from the zoom-in simu- 
lations with parameters from both WMAP1 (solid line) and WMAP3 (dashed 
line). The volume of the high-resolution zoom-in region is ~ 50 3 h~ 3 Mpc 3 . 

Sprin geT& Hernquisil d2003bl) : the wind efficiency 77 = 0.5, 
which measures the coefficient of the star formation that de- 
termines the mass outflow; the energy fraction from super- 
novae injected into the wind x = 0.25, wind free travel length 
L w = 20k pc, and a wind velo city V w = 41 8 km s" 1 . As demon- 
strated bv lCox et alj d2006al) . this wind model is able to repro- 
duce the starbursts as observed in Lyman Break Galaxies, and 
therefore is suitable to our study. 

We find that the impact of the starburst-driven wind on 
the quasar evolution is minor, as shown in Figure [14] Both 
the histories of star formation and black hole growth remain 
roughly the same as in the simulation without a starburst wind, 
only the amplitude is lowered by a factor of ~ 1.5. Similarly, 
the final masses of the black hole and the stars are reduced 
by roughly the same factor, but the quasar host is still on the 
MfiH-Mbuige correlation. The peak quasar phase is delayed 
to z ~ 6. Overall, the starburst wind affects the gas dynam- 
ics locally, but owing to the deep potential of the system, its 
impact on the process of qua sar formation is m inor. Our re- 
sults support the finding by ICox et all d2006al) that feedback 
from starburst-driven winds alone is ineffective at regulating 
the growth of the central black hole, so feedback from the 
black hole plays the dominant role in the formation and evo- 
lution of quasars. 

4.4. Abundance and Fate of Quasar Halos at z ~ 6 
Because we have so far simulated only one quasar in a vol- 
ume of l/T 3 Gpc 3 , we are not yet able to constrain the ex- 
pected abundance of quasars z ~ 6. As mentioned in § 2.3, 
at a given redshift, cosmological simulations with parameters 
from WMAP1 produce more massive halos than runs with 
WMAP3 owing to a larger value of erg. Figure Q3] shows the 
number of halos at z ~ 6 from the zoom-in runs with parame- 
ters from both WMAP1 and WMAP3. There are about three 
dozen halos with mass M > 1O 12 M in the WMAP1 run, 
while in the WMAP3 run there are only a handful of such ha- 



los. However, since in our picture the quasar activity depends 
not only on the halo mass, but also on the merging history, 
an accurate estimate of the quasar abundance and luminos- 
ity function would require hydrodynamical simulations of all 
the quasar candidates in a large box, which are currently un- 
available. Nevertheless, all conditions being equal, the change 
from the WMAP3 parameters would produce fewer luminous 
quasars at z ~ 6. This suggests that in a WMAP3 cosmology, 
the quasar observed with the largest redshift, Jl 148+5251, 
might have formed in a slightly higher overdensity peak than 
that we have presented here. In that event, if the WMAP3 de- 
termination of a% were correct, we would need to identify a 
rarer density fluctuation to match J 1 148+5251 at its observed 
redshift. However, this does not change our conclusion that 
the most distant and luminous quasars can form from hierar- 
chical galaxy mergers in the ACDM cosmology. 

Imaging surveys of J 1 148+5251 show that the re is no other 
luminous quasar from the same epoch in the fi eld dCarilli et al.l 
I2004t IWhite et alJl2005t IWillott et al.ll2005l) . In our simula- 
tions, around the peak of quasar activity at z ~ 6.5, there are 
no other halos of mass > 10 12 M Q within a few Mpc of this 
quasar. However, the numerous major mergers this halo ex- 
perienced prior to the peak quasar activity demonstrates that 
this region was once highly clustered with massive halos, but 
they merged to become a bigger one by z = 6.5. 

As seen from FigureO this quasar halo will undergo a hand- 
ful of major mergers at a later time from z ~ 4-1, and eventu- 
ally end up as a cD-like galaxy at the center of a rich cluster. 
Since we do not follow hydrodynamically the evolution of the 
quasar at z < 4, the physical conditions of these mergers re- 
main undetermined. It is not clear whether this halo would 
experience more episodes of starburst or quasar activity later 
on during these mergers. Therefore, the final black hole mass 
and other properties of this quasar at the present day are de- 
ferred to future simulations that follow its evolution to z = 0. 

4.5. Galaxies in the Epoch of Reionization 

The epoch of reionization (EoR) is an important landmark 
event in cosmic histo ry that constrains the f ormation of the 
first luminous objects dLoeb & Barkana 2001). The recent re- 
sults of WMAP 3 indicate that the Universe was 50% reion- 
ized at z ~ 9.3 (Pag e et al.|[2006t ISpergel et ail 120061) while 
studies of Gunn-Peterson absorption (|Gunn & Peterso n! 1965b 
suggest th at reionization began as early as z ~ 14 and ended 
at z ~ 6 dFan et al .1 12006). At present, it is believed that the 
reionization sources are star-forming galaxies since there are 
insufficient quasars at z > 6 as indic ated by the steep quasar 
luminosity function dFan et al.ll2006l) . 

The galaxy progenitors of the quasar in our simulations 
underwent extreme and prolonged starbursts before z ~ 
6.5. Less extreme galaxies in this epoch may also have 
vigorous star formation histories. Detecting these galax- 
ies and determining their contribution to reionization will 
be crucial to understanding the EoR (iHernquist & Sp ringel 
| 2003b iBarton et al.ll2004t iDave et al.ll2006h . As reviewed by 
iHu & Cowi e (2006), recent observations u sing both broad- 
band colors (e.g., [Dickinso n et al.l l2004t lYan et all 120051; 
iBunker et all l2004tnBouwens et all [20041: iGiavalisco et al l 
1 20041; [Egami et al.l 120051; lEvles et al.l 120051; iMobasher etall 
120051; lYan et all 120061; lEvles et al l 120061) and narrow-band 
Lyq emission (e. g. jHu et al .1120021: iMalhotra & Rhoadsll2004t 
IStern et aI1l2005l) h ave detected ~ 500 g alaxies at z ~ 6 and 
a handful at z > 7 dBouwens et al.ll2005l) . The low luminos- 
ity density of galaxies currently detected at z > 7 seems in- 



22 



Li et al. 



sufficient to reionize the Universe. However, ongoing surveys 
with the HST and Spitzer telescopes, and f uture missions such 
as the Dark Ages z Lya Explorer (DAzLE. lHorton et al.l2004l) 
and th e James Webb Space Telescope (JWST, Win dhorst et al.l 
20061) will search deeper and further for high-redshift objects, 
and may eventually unveil ionizing sources in the EoR. 

5. SUMMARY 

We have presented a model that accounts for the SMBH 
growth, quasar activity and host galaxy properties of the most 
distant quasar observed at z = 6.42, by following the hierar- 
chical assembly of the quasar halo in the standard ACDM 
cosmology. We employ a set of multi-scale simulations that 
include large-scale cosmological N-body calculations and hy- 
drodynamic simulations of galaxy mergers, and a recipe for 
black hole growth self-regulated by feedback. We first per- 
form a coarse N-body simulation in a volume of 1 h~ 3 Gpc 3 to 
identify the largest halo at z = 0, which is assumed to be the 
descendant of the earliest luminous quasar. We then "zoom 
in" on the halo and resimulate the region with higher reso- 
lution sufficient to extract its merging history starting from 
very high redshift. The largest halo at z ~ 6 reaches a mass 
of ~5.4x 10 12 /i _1 Mq through 7 major mergers betweenz~ 
14.4-6.5. These major mergers are again re-simulated hydro- 
dynamically using galaxy models and recipes for star forma- 
tion, SMBH growth, and feedback. 

We find that the quasar host galaxy builds up rapidly 
through gas-rich major mergers, with star formation rates 
up to 10 M Q yr _I , reaching a stellar mass of - 10 12 M Q at 
z ~ 6.5. The black holes grow through gas accretion under the 
Eddington limit in a self -regulated manner owing to their own 
feedback. As the galaxies merge, the black holes coalescence 
to form a dominant black hole, reaching a peak accretion rate 
of ~ 20M© yr" 1 and a mass of M B h ~2x 10 9 M q at z ~ 6.5. 
Feedback from black hole accretion clears away the obscuring 
gas from the central regions, making the quasar optically vis- 
ible from z ~ 7.5-6. At the peak of the quasar phase, the star 
formation rate, metallicity, black hole mass, as well as quasar 
luminosities of the simulated system are consistent with ob- 
servations of Jl 148+525 1 . 

Our results demonstrate that rare and luminous quasars at 



high redshifts can form in the standard ACDM cosmology 
through hierarchical, gas-rich mergers, within the available 
cosmic time up to the early epoch of z ~ 6.5, without requiring 
exotic processes. Our model should also provide a viable for- 
mation mechanism for other distant, luminous quasars. More- 
over, we predict that quasar hosts at high redshifts follow sim- 
ilar Mbh - Mbuige correlation observed locally as a result of 
the coeval evolution of the SMBHs and host galaxies. Bet- 
ter measurements of black hole masses and host properties 
with future observations will therefore be crucial to test our 
prediction. Furthermore, we predict that the progenitors of 
the distant quasars undergo strong and prolonged starbursts 
with rates ~ 10 3 M Q yr _1 at higher redshifts z > 8, that would 
contribute to the reionization of the Universe. Detecting these 
early galaxies and unveiling the epoch of reionization will be 
an important goal of current and future missions in observa- 
tional cosmology. 
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